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PREFACE 


Three-dimensional fully self-consistent computer models were used to deter- 
mine the evolution of galaxies consisting of 100 000 simulation stars. One 
series of computer experiments used initially balanced flattened galaxies with 
an exponential radial density variation. Comparison of two-dimensional (infin- 
itesimally thin disk) simulations with three-dimensional (disk with finite 
thickness) simulations showed only a very slight stabilizing effect due to the 
additional degree of freedom. The addition of a fully self-consistent, nonro- 
tating, exponential core/halo component resulted in considerable stabilization. 
The most pronounced instabilities present were those due to the Jeans' insta- 
bility in the outer regions of the disk, while at the same time a relatively 
slowly growing bar instability appeared. A second series of computer experi- 
ments was performed to determine the collapse and relaxation of initially spher- 
ical, uniform density and uniform velocity dispersion stellar systems. The 
evolution of the system was followed for various amounts of angular momentum 
in solid body rotation. For initially low values of the angular momentum satis- 
fying the Ostriker -Peebles stability criterion, the systems quickly relax to 
an axisymmetric shape and resemble elliptical galaxies in appearance. The max- 
imum flattening for these systems is equivalent only to an E2 system. For 
larger values of the initial angular momentum, bars develop and the systems 
undergo a much more drastic evolution. The apparent rotational and random 
velocities of the bar systems are very sensitive to the viewing direction. An 
additional complication is the frequent misalignment of the apparent major axis 
with the direction that reflects the maximum rotation. 
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INTRODUCTION 


In recent years, large-scale N-body computer simulations (refs. 1 and 2) 
have become an important tool in investigating the structure of spiral galax- 
ies, especially in determining the development of large-scale instabilities 
resulting in spiral and bar formation. Until recently, most of these simula- 
tions used essentially two-dimensional models with the "stars" confined to the 
plane of the galactic disk (refs. 3 and 4) . These simulations have shown that 
the disks of stars tend to develop fast-growing nonaxisymmetric instabilities 
resulting in bar formation. The bar instabilities occur even for velocity dis- 
persions that are considerably larger than those found in the solar neighbor- 
hood or those predicted by Toomre (ref. 5) as being locally stabilizing. Global 
stability studies of disks of stars have been primarily numerical large-scale 
N-body simulations. Some limited analytical work has been done primarily for 
uniformly rotating disks (refs. 6 and 7) , but generally linear stability anal- 
yses were used in these analytical studies of disks of stars. 

Any spiral structure in computer-generated galaxies is generally short 
lived and the final state is a rotating bar. The bar thus obtained rotates 
more slowly than the stars. For a typical case investigated by Hohl (ref. 4) , 
the bar rotates at 2.25T and the stars rotate at 1.5T where t is the rota- 
tional period of the initial disk. 

It has been argued by Ostriker and Peebles (ref. 8) that core/halo compo- 
nents have a stabilizing effect on galaxies and result in longer lived spiral 
structures. However, numerical experiments with large fixed stellar components 
representing the core/halo component (refs. 9 and 10) show that a multiarmed 
spiral structure develops and persists for many rotations but only in an evolv- 
ing manner. That is, the spiral structure is either wound up into a tight pat- 
tern. or it is wound up and then reappears again. A recent study of the effect 
of fixed core/halo components using two-dimensional N-body simulations (ref. 11) 
does show that the bar instability is indeed inhibited by a sufficiently large 
fixed component. 

The purpose of the present study of initially balanced flattened galaxies 
is twofold. First, we want to determine the effect of a self-consistent (rather 
than a fixed nonresponding) core/halo component. This will show whether there 
are any instabilities (such as "two-stream") or other important interactions 
present that may be suppressed with a fixed nonresponding core. Second, we 
want to determine the effects of finite thickness of the disk and of three- 
dimensional essentially spherical core/halo components. 

The second part of the present study, using initially spherical collapsing 
galaxies, is aimed at resolving the dynamics of the formation of galaxies. It 
is likely that galaxies arise directly from the gravitational collapse and sub- 
sequent relaxation of portions of the expanding universe. After initially col- 
lapsing, the system expands and undergoes additional oscillations with ampli- 
tudes and damping depending on the dissipative mechanisms operating. According 



to this script, elliptical galaxies result if the gaseous component is largely 
exhausted by the end of the first collapse. The various physical processes 
involved in this simple collapse picture, as well as some nonconventional alter- 
natives, have been discussed in the reviews by Larson, by Jones, and by Gott 
(refs. 12 to 14). For situations in which very little gas is left by the time 
of first collapse, a collisionless system of stars appears to be an appropriate 
model. Numerical experiments studying the effect of dynamical mixing in such 
systems with spherical symmetry were carried out by Bouvier and Janin (ref. 15) 
and by Henon (ref. 16). 

Less restricted simulations aimed specifically at modeling the formation 
of elliptical galaxies by collisionless stellar dynamics were performed by Gott 
(refs. 17 and 18), who included rotation and relaxed the imposed symmetry condi- 
tions from spherical to axial. These uniformly rotating models produced flat- 
tened, oblate systems. In later work, Gott (ref. 18) added the effects of cos- 
mological infall and tidal interactions and obtained good agreement with the 
ellipticity profile and light distributions of ellipticals. However, Gott's 
calculations have remained somewhat suspect since his method did not permit the 
development of the fierce bar instabilities and the associated large-scale 
momentum transfer that have been so common in simulations of disk galaxies 
(refs. 4, 8, and 11). 

Recently Miller (ref. 19) and Miller and Smith (ref. 20) have produced 
some fully three-dimensional simulations of the collapse of stellar systems. 
Their earlier examples produced prolate bars rather than oblate systems like 
those obtained by Gott. This striking difference, however, seems attributable 
more to the colder initial conditions used by Miller and by Miller and Smith 
than to their removal of the constraint of axial symmetry. 

But even leaving aside this concern about possible suppressed nonaxisym- 
metric instabilities, the merits of the oblate models of ellipticals have been 
called into question by the recent observations of the rotational and random 
velocities of elliptical systems by Illingworth (ref. 21) and by Peterson 
(ref. 22) . As discussed by Illingworth, data for numerous elliptical galaxies 
seem to indicate that the ratios of rotational to random velocities are typi- 
cally a factor of 2 or 3 smaller than those predicted by the oblate models. 

It may well be that more exotic models of the formation of elliptical gal- 
axies are required. Nevertheless, it seems a pity to give up too soon on the 
simple collapse model, especially since the range of initial conditions covered 
by the fully three-dimensional calculations of Miller and of Miller and Smith 
is limited. We report here some extensive results of three-dimensional models 
of the collisionless collapse and relaxation of rotating stellar systems with 
initial random velocities patterned after those of Gott (ref. 17). This per- 
mits an assessment of the importance of his requirement of axial symmetry. At 
the same time, these initial conditions cover a markedly different range from 
those of other three-dimensional experiments. 
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SYMBOLS 


C 

c 

G 

H 

h 

I 

K 

Mo 

m 

N 

n 

P 

Q 

R 

r,4>,0 

T 

T cir 

T rand 

t 

U 

V 

V m 


constant, 3M g /(4TTr 2 P(0,0)) 
velocity, defined in equation (5) 
gravitational constant 
Green's function 

defines z-dimension of n x n x h active array 

moment of inertia 

gravitational field 

mass of galaxy 

mass of Sun 

mass of simulation star 

dimension of array used in potential calculation 

defines n x n x h active array, N/2 

angular momentum 

= a r /tf r ,mi n 

radius of galaxy 

spherical coordinates 

random kinetic energy 

kinetic energy in rotation 

random kinetic energy 

time 

■ V m/<M°) 

velocity 

peak rotational velocity 
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w 


x,y,z 


e 

? 

K 


y 

p 

P (0) 

a 

°r ,min 
a v 

4 > 

a 


total potential energy 
Cartesian coordinates 
ellipticity 

radius, x 2 + y 2 + (7z/5 ) 2 

epicylic frequency 

surface mass density 

volume mass density 

mass density at center of galaxy 

velocity dispersion 

velocity dispersion defined by equation (3) 
line-of-sight velocity dispersion 
gravitational potential 
angular velocity of spherical model 
- \/GMg/R 3 


w angular velocity of disk model 

(jJ 0 angular velocity of cold balanced disk 

Subscripts: 

i,j,k summation indices 

r,4>,9 radial, azimuthal, and colatitudinal component 

x,y,z x-, y-, and z-component 

summation indices 


Notation: 

Fourier transformed quantity 


MODEL 

The model used for the present galaxy simulations consists of 100 000 rep- 
resentative stars of equal mass that move inside an active 64 x 64 x 16 (or 
33 x 33 x 33 ) three-dimensional array of cells. For the disk simulations the 
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stars are confined to the plane of the disk represented by a 64 * 64 active 
array. The number of stars inside a particular cell defines the mass density 
at the center of each cell. Fast Fourier transform methods are used to obtain 
the gravitational potential at the center of each cell for a given mass density 
distribution. Thus the potential calculation actually requires an array larger 
than the 33 * 33 x 33 active array so that the periodicity associated with the 
Fourier transform can be eliminated. The force acting on a particular star is 
determined by trilinear interpolation from the gravitational field at the sur- 
rounding eight cell centers. After the force acting on a star is determined, 
its position and velocity are advanced by a small time step using the time- 
centered leap-frog method. Details of the two-dimensional model are given by 
Hohl and Hockney (ref. 2) and by Hohl (ref. 23). The extension of the model 
to three dimensions is described in the appendix. Collisional effects in this 
model become important only after hundreds of rotations (ref. 24) . 


DISK-CORE SIMULATIONS 

Observational evidence (refs. 25 to 27) indicates that the luminosity (and 
presumably the density) in the outer regions of many spiral and SO galaxies 
decreases exponentially with radius. Also, previous simulations (ref. 4) showed 
that initially unstable stellar disks evolved into stable systems with radial 
density variations that closely approximated the sum of two exponentials. The 
inner exponential with a scale length of about 1 kpc describes the nonrotating 
or slowly rotating spheroidal, or core, component and the remaining exponential 
with a scale length of about 8 kpc describes the extended disk population. 

Thus, it seems reasonable to use an exponential density variation for the disk 
in the present computer simulations. Similarly, the central core used is 
described by an exponential density variation. 

Figure 1 illustrates the evolution of a disk of 100 000 stars with an ini- 
tially exponential surface mass density distribution P ( r ) = y ( 0 ) e _r / 2 with 
a cutoff at r = 10 kpc. The initial angular velocity of the disk was obtained 
from 


u)2 = Uq 2 + _J — [y(r) o r 2 (r) ] 

r y (r) 3r 



with 


°<J> (r) 


Mr) 

2 o) Q (r) 


a r (r) 


( 2 ) 


Here, ( 0 o (r) is the angular velocity required to balance the cold (zero veloc- 
ity dispersion) disk, co(r) is the actual angular velocity, and Mr) is the 
epicyclic frequency. The initial value of the radial velocity dispersion o r 
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was taken to be that determined by Toomre (ref. 5) as the minimum required to 
stabilize all local axisymmetric instabilities. 


a r (r) = a r>min = 3.36G u(r)/K(r) 


( 3 ) 


The time t is given in rotational periods (2n/(D 0 ) of the cold disk at a radius 
of 5 kpc, that is, halfway to the edge of the initial disk. 

As expected (refs. 4 and 9) , only the small-scale instabilities are pre- 
vented by a r = CT r/in i n and the system quickly forms a two-arm spiral which 
eventually tends to evolve into a rotating bar. The evolution of the azimuthally 
averaged radial density variation for this system is shown in figure 2. As pre- 
viously observed (refs. 4 and 9) , the eventual density variation approaches one 
which can be closely approximated by the sum of two exponentials. One exponen- 
tial describes the central core component and the other describes the extended 
disk. The evolution of the radial and azimuthal velocity distributions as a 
function of radius is shown in figure 3. It is obvious that the system has not 
reached an equilibrium condition. Note that there are still large outward veloc- 
ities in the outer regions of the system. This can be seen more clearly in fig- 
ure 4 where the evolution of the azimuthally averaged, mean radial velocity is 
plotted as a function of radius. This figure displays considerable systematic 
radial flows. The radial and azimuthal velocity dispersions displayed in fig- 
ure 5 show that there is some heating near the center and that the velocity dis- 
persion for stars expanding into the extended disk component increases consider- 
ably. The increases in the central mass density and heating near the center 
result in large changes in the rotation curve. Figure 6 shows the evolution 
of the rotation curve. 


v <{> 0 = r w o( r ) = >/rK j 


and the corresponding mean rotation of the stars <V(jj>. In the central region 
of the system the galaxy is now balanced to a large extent by increased pressure 
due to the higher velocity dispersion. The mean rotation has been reduced in 
this region by transferring momentum to larger radii. This is illustrated in 
figure 7 which shows that as the bar instability develops (compare with fig. 1), 
much of the angular momentum is transferred to larger radii. 

Various other diagnostics have been performed on the system. For example, 
figure 8 shows the evolution of the moment of inertia I and the angular momen- 
tum P divided by their values at t = 0. As can be seen, P is conserved in 
the simulation, but I is still increasing at a nearly linear rate after three 
rotations. The evolution of various components of the total kinetic energy 
divided by the total potential energy is shown in figure 9. The components 
T r and T* represent the random kinetic energies due to the velocity disper- 
sions in tne r- and (^-directions, respectively, while T c i r is the kinetic 
energy of rotation. Note that the ratio of the kinetic energy in rotation to 
the absolute value of the total gravitational energy of the system is approach- 
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ing the value 0.14 predicted by Ostriker and Peebles (ref. 8) for stability. 
This indicates that considerable heating of the system occurs. 

One of the aims of the present study is to determine the effect of adding 
the third degree of freedom by allowing the exponential disk to have a finite 
thickness. Using again an exponential projected surface mass density variation 
y(r) = y(0) e ~ r// the stars are now distributed in the z-direction according 
to the one-dimensional distribution sech^ (z/C) . The parameter C is deter- 
mined from p (0 ) = 3Mg/4nR^C (ref. 28) where Mg is the total mass of the 
galaxy and R = 10 kpc is the radius of the disk. The central thickness of 
the disk is 2 kpc and the density is cut off at z-| given by 



sech^ 



0.1 


(4) 


The radial and azimuthal velocity components are determined in a manner similar 
to that for the infinitesimally thin disk and the z-component of the velocity 
dispersion is determined by a force balance in the z-direction. Note also that 
all initial velocities are truncated so that stars have kinetic energies no 
greater than those which would allow them to reach the boundary of the system 
in the gravitational potential at t = 0. 

Figure 10 shows two side views of the initial disk and the evolution for 
up to 3 rotations. The differences in the projections as the evolution proceeds 
are due to the bar formation. Note also the rapid expansion in the plane of the 
disk. This is the result of the bar instability, as shown in figure 11 which 
gives the evolution of the disk projected in the x-y plane. Note that the evo- 
lution is very similar to that shown in figure 1 for the thin disk. The simi- 
larities in the evolution of the two systems are even more accentuated by com- 
paring the evolution of the surface mass density distribution and the velocity 
distribution in figures 12 and 13 with those in figure 2 and 3. For complete- 
ness figures 14 to 19 show the evolution of the various other parameters for 
comparison with figures 4 to 9. For the finite thickness system, the additional 
variable of z-component velocity or kinetic energy remains small compared with 
the r- or ^-components. Thus, one would expect little difference in the evolu- 
tion of the finite thickness disk when compared with the infinitesimally thin 
disk. 


As shown in figures 1 and 11, exponential disks with velocity dispersion 
(Q « 1) are violently unstable to the bar-forming node. Previous work (ref. 11) 
with a superimposed fixed (non-self-consistent) central mass distribution indi- 
cated a stabilizing effect toward the bar-forming instability. A more realistic 
simulation is to allow core-disk interaction; thus, we are presently interested 
in the stabilizing effects of a completely self-consistent core, or spheroid, 
component. Again, the effect is investigated for both the infinitesimally thin 
disk (two-dimensional) and for the three-dimensional disk. 


For the core-disk system, 50 percent of the mass (50 000 stars) is con- 
tained in the nonrotating core and the remaining mass (50 000 stars) is con- 
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tained in the disk. The disk component is again given .the surface density vari- 
ation Pdisk( r ) = ^disk(°) e -r /2, whereas the initial nonrotating core compo- 
nent is given a density variation Ucore( r ) = Mcore(O) e"^ r . Note that the 
disk and core densities are cut off at r = 10 kpc and r = 3.5 kpc, respec- 
tively. The initial velocity dispersion and rotation of the disk are obtained 
by again using equations (1), (2), and (3) with y = Mdisk* Similarly, as 
before, the z-dimensions of the disk are determined from equation (4) . The 
initial velocity dispersion of the nonrotating core component was obtained by 
taking = a r and simply balancing the core in the presence of the disk. In 
order to ensure that the core component was in a stable equilibrium state at 
the start of the core-disk simulation, the core was allowed to evolve for sev- 
eral rotational periods (2 tt/o) 0 at 5 kpc) with the disk stars held fixed. Start- 
ing from these initial conditions, the system evolved as shown in figure 20. 

Note that even though a two-arm spiral structure still forms, the system as a 
a whole evolves in a much less violent manner than that displayed in figure 1 . 
This can also be seen in figure 21 which shows the evolution of the surface mass 
density for both the core and the disk components. Nevertheless, the spiral 
instability in figure 20 is leading toward the formation of a bar structure. 

Note also that with the exception of a slow outward diffusion of stars near the 
edge, the core remains essentially stationary, while the disk component displays 
the outward shift of mass generally associated with bar formation. Figure 22 
shows the evolution of the radial and azimuthal velocity distributions as a 
function of radius, when comparing the results with figure 3, two features 
stand out. One feature is, of course, the large-velocity stars in the core near 
the center and the other is the much more violent instability at larger radii in 
figure 3. The mean radial velocity shown in figure 23 also displays a consider- 
able reduction in the net radial flow when compared with figure 4. Similar 
information is contained in figure 24 which displays the evolution of the radial 
and azimuthal velocity dispersions for the core and disk components. Note the 
sharp increase in the radial velocity dispersion at r » 3 kpc which is associ- 
ated with a marked reduction in the angular momentum of the disk in this region. 
This outward shift of the angular momentum is displayed in figure 25 which shows 
the evolution of the angular momentum distribution. In general, the simulations 
show that the formation of bars or two-armed spirals results in moving angular 
momentum outward to larger radii. The evolution of the rotation curve displayed 
in figure 26 shows that for up to 3 rotations the disk component continues to 
have a high rotational velocity in the central region of the system. 

As shown in figure 27, the addition of the core component markedly reduced 
the rate of increase in the moment of inertia. Similarly, the ratios of kinetic 
energy to potential energy shown in figure 28 change considerably less than 
those of the disk without a central core component. This reflects the sizable 
differences in the initial values of these ratios. 

The final system investigated in this series is a three-dimensional expo- 
nential disk with a three-dimensional core, or spheroid, component. The spa- 
tial distribution of the stars for the disk component is obtained as was done 
for the disk shown in figures 10 and 11 except that now the disk contains only 
50 000 stars. For the nonrotating central core, the density is given by 
P (C) = P (0) e“2?, where Xp- = x^ + y2 + (7z/5)2. The density is cut off at 
£ = 7. Thus, the central core or spheroid has an axis ratio of 7:5. Again 
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the Gaussian velocity dispersion for the core is obtained by a simple balance 
of the self-gravity of the total system. The velocity dispersion for the disk 
component is generated as was done for the system shown in figure 11. Before 
initiating the simulation of the combined core-disk system, the core was allowed 
to evolve for several rotations (with fixed disk stars) to ensure that no insta- 
bilities or other problems associated with the core component were present. 

Figure 29 shows the evolution of the system perpendicular to the equatorial 
plane. Note the remarkable stability of the system when compared with the disk 
without the central core in figure 10. Figure 29 does show a slight bulging of 
the disk at t = 3.0 indicative of spiral-type instabilities. The evolution of 
the system in the equatorial plane is shown in figure 30 and displays the devel- 
opment of a comparatively weak spiral structure initiated by Jeans' type insta- 
bilities in the outer region of the disk. Although it is not apparent from fig- 
ure 30, an analysis of the azimuthal density variation for the disk stars at 
t = 3.0 shows that a bar structure is forming. For example, in the region from 
2 to 5 kpc, which is near the center away from the visible spiral structure dis- 
played in figure 30, the peak density along the bar is 2.5 times the density 90° 
away from the major bar axis. It should be noted that because of the allowed 
initial relaxation, the core components of the two core-disk systems investi- 
gated here are expected to closely satisfy the collisionless Boltzmann equation. 
The same is not necessarily true for the disk component since satisfying equa- 
tion (1) only ensures a balance of forces at t = 0. Also, for a stellar disk, 
a r = ^r, min <3°e s not ensure stabilization of global nonaxisymmetr ic instabili- 
ties (refs. 29 and 30) . However, since one would hardly expect nature to gener- 
ate a galaxy initially in an exact stable stationary state and since we are 
interested in the future development of instabilities and the final state toward 
which the system evolves, an exact stationary and stable initial state is not 
necessary. 

The evolution of the radial, azimuthal, and axial velocity distributions 
as a function of radius is shown in figure 31 . Note that now there is even 
less change visible than for the two-dimensional system shown in figure 22. 

This is illustrated more clearly in figure 32 which shows very little net mass 
motion in the radial direction. 

The evolution of the azimuthally averaged projected surface mass density 
for the three-dimensional disk-core system is shown in figure 33 and is nearly 
identical to that of the two-dimensional disk-core system shown in figure 21. 
Note that there is very little change in the density for the core with the 
exception of a slight outward diffusion near the edge. Azimuthally averaged 
values of the total density variation in the z-direction are shown in figure 34 
for various values of r. Note that the density follows an essentially exponen- 
tial variation. Sane of the fluctuations shown may be due to the relatively 
small sampling volume used. The evolution of the radial velocity dispersion 
shown in figure 35 indicates that (as expected) the velocity dispersion for the 
two-dimensional core (shown in fig. 24) is higher. Also, the large increase in 
the velocity dispersion of the thin disk near r = 3 kpc does not occur for 
the three-dimensional disk. Associated with this is the fact that there is very 
little change in the radial angular momentum distribution (fig. 36) during the 
evolution of the three-dimensional disk-core system, whereas considerable out- 
ward shift of angular momentum occurs for the two-dimensional disk-core system 
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(fig. 25) . These results indicated that the global bar instability is much 
weaker for the three-dimensional system than for the two-dimensional system, 
as is obvious by comparing figures 20 and 30. 

The evolution of the moment of inertia and of the kinetic energy ratios 
for the three-dimensional disk-core system is shown in figures 37 and 38. As 
can be seen, there is little change in the value of the various components dur- 
ing the evolution. Note that the ratio of the kinetic energy in rotation to 
the total potential energy of the system is slightly higher than the value of 
the 0.14 predicted for stability by Ostriker and Peebles (ref. 8). Also, the 
moment of inertia increases by only about one-third of that shown in figure 27 
for the two-dimensional system. As was the case for all four systems investi- 
gated, the angular momentum was conserved to a high degree of accuracy. 


SPHERICAL COLLAPSE SIMULATIONS 

The initial locations of the stars for this series of simulations were 
obtained by distributing them at random inside a sphere of radius R = 10 kpc 
to produce a uniform mass density. The total mass of the system was taken to 
be Mg = 2 x IO^Mq. Solid body rotation about the vertical axis with angular 
velocity 


a 


o 



29.3 km- sec -1 -kpc' 1 


balances the cold (zero velocity dispersion) system in the equatorial plane. 
The initial velocities of the stars were selected from an isotropic Maxwellian 
distribution superimposed upon solid body rotation with angular velocity 

= ojSIq. The amounts of rotation and the radial velocity dispersion a r used 
for each of the five models are listed in table I. Also given are the initial 


TABLE I.- INITIAL CONDITIONS 


- 

Model 

I 

II 

III 

IV 

V 

n/ft 0 

o r , km/sec . . . 

T r and/ 1^1 • • • 

T cir /|w| . . . 

0 

92.8 

1/4 

0 

0.500 

92.8 

1/4 

1/12 

0.707 

92.8 

1/4 

1/6 

0.866 

65.7 

1/8 

1/4 

1 .159 
41 .5 
1/20 
4/9 
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ratios of the random kinetic energy T ran ^ and the rotational kinetic energy 
T cir to total gravitational energy W. The random kinetic energy is given 
by 


10 5 

T rand a ^ ( c r,i + c 9,i + c <(),i) 
i=l 

where m is the mass of a star and 

c? = (V r - <V r » 2 ' 

C0 = (Vq - <v 0 » 2 > 
c| = (V(j, - <v 4 ,» 2 > 

Note that initially (at t = 0) 

< v r> ' < v 0> = 0 < v 4>> = ( r sin 

The rotational kinetic energy is given by 


(5) 


T 


cir 




< v 0,i> 


The first four models in table I have precisely the initial conditions used by 
Gott (ref. 17). 


Main Features 

The collapse and relaxation of the five models are illustrated in fig- 
ures 39 through 43. For reasons of computational efficiency, the snapshots 
of the three orthogonal projections shown at each time in these figures were 
actually taken a time step apart. This accounts for the apparent lack of sym- 
metry at certain stages. Since the initial configurations of each of these 
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models were far from equilibrium, the initial evolution is characterized by the 
collisionless violent relaxation process described by Lynden-Bell (ref. 31). 

In this process the distribution function tends, albeit incompletely, toward 
a Maxwellian form through the influences of the rapidly fluctuating gravita- 
tional field. The effectiveness of this process is expected to decrease with 
radial distance from the center since the duration of the rapidly fluctuating 
field is a decreasing fraction of the orbital time scales. 

For models I, II, and III (the hotter and less rapidly rotating ones), this 
violent relaxation is the principal effect which determines the final, equilib- 
rium configuration. On the other hand, in models IV and V an additional effect 
occurs, the development of a pronounced bar-shaped structure. For this reason, 
model IV is selected for special attention; not only is its evolution followed 
through 5 rotations rather than 3, but also the full time development of its 
properties is displayed rather than just the final properties, as for the 
remaining models. Because of difficulties in preserving the resolution of the 
nonaxisymmetric structure during the photographic processing, models IV and V 
in figures 42 and 43 are actually shown twice. Figure 42(a) shows model IV in 
three projections of 100 000 stars, each. For figure 42(b), the projections are 
printed much lighter and the x-y projection contains only 50 000 stars with 
100 000 stars for the x-z and y-z projections. The same procedure is used for 
the two displays of model V in figures 43(a) and 43(b). 

From the final frames of figures 42 and 43, it is apparent that a bar sys- 
tem appears quite flattened for certain viewing angles. The degree of flatten- 
ing can be estimated from projected isodensity contours (see the discussion of 
fig. 66 subsequently) . The eccentricities of the contours in the bar models 
decrease with increasing radius. The axial ratios of 1.0 to 0.39 to 0.54 and 
1.0 to 0.19 to 0.28 for models IV and V were obtained from the innermost iso- 
density contours. The vertical axis is the last item listed for each model. 
Thus, even an apparent E8 system can be obtained from a collisionless model 
without invoking any dissipation. Notice also that the shortest of the three 
axes for the bar models lies in the equatorial plane rather than along the axis 
of rotation. These are thin bars rather than flat ones, as indeed is evident 
from figures 42 and 43. On the other hand, the isodensity contours for the 
axisymmetric models II and III are noticeably more eccentric in the outer 
regions than in the center. The maximum eccentricities found are 0.18 for 
model II and 0.25 for model III. 

Figures 44 through 47 are useful in describing the general features of the 
collapse of the models. Both the moment of inertia and the various kinetic 
energy ratios reflect the oscillations occurring during the collapse phase. 
Although the moment of inertia oscillates only near its initial value for 
model I, for models II to IV it displays a steady increase after undergoing one 
or more oscillations reflecting the initial collapse in the horizontal direc- 
tion, that is, in the direction perpendicular to the rotation axis. Model V, 
of course, initially expands in the horizontal direction since the centrifugal 
forces alone are sufficient to overcome the gravitational attraction. The 
eventual steady increase in the moment of inertia is due to a comparatively 
small number of escaping stars rather than to any failure to achieve equilib- 
rium. The constancy of the total angular momentum, which is shown in figures 44 
and 45, is reassuring evidence of the accuracy of the numerical simulation. 
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For models lv and V (shown in figs. 42 and 43) , the first indications of 
the developing bar occur at t = 1 .5 and t = 0.75, respectively, in the form 
of vertical "flares" at both ends of the two side projections. In subsequent 
frames the nonaxisymmetric structure is clearest in the equatorial (top) pro- 
jection. With the exception of the vertical flares, the transitory development 
of the bar is similar to the bar formation process in the three-dimensional 
simulations of equilibrium, but unstable, flattened galaxies discussed earlier. 
Bar formation appears to become significant during or after the second collapse. 
For example, in figure 43, the peak of the second collapse, as indicated by the 
maximum of the z-component of the kinetic energy, occurs shortly before 
t = 0.75. At that time, large-scale axisymmetric structures have appeared. 

Prior to the emergence of this instability it is evident that the evolution is 
dominated by simple, axially symmetric collapse dynamics. 

The amount of time taken by the initial collapse may be compared with that 
required for the collapse of the corresponding pressure-free, uniformly rotating 
Maclaurin spheroid. This simple model was described by Lynden-Bell (ref. 32) 
and has been used for comparison with stellar dynamic results by Gott and Thuan 
(ref. 33) and by Miller and Smith (ref. 20). It predicts that the end of the 
first collapse in the vertical direction is reached at t = 0.18, 0.19, 0.20, 
0.20, and 0.21 for models I, II, III, IV, and V, respectively. The location of 
the first peak in the vertical kinetic energy ratios (figs. 46 and 47) suggests 
that the initial collapse is completed near t = 0.25 for models I to III and 
near t = 0.20 for models IV and V. Notice that in the four rotating models 
the horizontal components of the kinetic energy reach their first maximum dis- 
tinctly later than the vertical components, as expected. Di Fazio and 
Occhionero (ref. 34 ) add pressure and relaxation effects to the Maclaurin 
spheroid model. Both of these effects increase the estimated collapse time, a 
result which appears needed for the three hottest models. However, a detailed 
comparison seems unwarranted in view of the uncertain relation between stellar 
velocity dispersion and the specific form of pressure employed by them. 

The kinetic energy ratios suggest that models I to III have fairly well 
relaxed by t = 1 . They also suggest that the collapse phase of model IV is 
largely completed by then. There appears to be a distinct interval between the 
end of the collapse process and the onset of the bar formation. For model V, 
however, it is more difficult to tell when the collapse stops and the bar forma- 
tion begins. The duration of the collapse in models I to IV conforms to Gott's 
result showing that the violent relaxation occurs in roughly 3 collapse times. 

In the present units, the reference collapse times used by Gott are 0.27, 0.32, 
0.40, and 0.36 for models I to IV and 0.49 for model V. A further comparison is 
available of the total maximum kinetic energy achieved during the collapse: 
for models I to IV, the maximum values of T/|w| are 0.65, 0.61, 0.57, and 0.58 
compared with 0.66, 0.62, 0.58, and 0.60 obtained in the axially symmetric cal- 
culation by Gott. 

A more detailed critique of the axially symmetric simulations can be made 
by comparing the distribution of the stars with respect to their energy. The 
results for the three-dimensional simulations are shown in figures 48 and 49. 

The final distributions of non-bar-forming models I to III are essentially the 
same as those obtained previously by Gott. The basic feature of a separation 
into a low energy core and a high energy halo is apparent. Even the locations 
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and strengths of these two components correspond closely to those found by Gott. 
As might be expected from the rapid equilibration of the kinetic energy ratios, 
the energy distributions of the first three models at t = 3 differ only mar- 
ginally from those at t = 2. Although the similarity of the present energy 
distributions to those obtained by Gott does not preclude a redistribution of 
angular momentum, there is no evidence that any redistribution occurs. Indeed, 
for each of these three axisymmetric models, the initial and final angular 
momentum distributions are virtually indistinguishable. 

The occurrence of a bar in models IV and V has a marked effect upon the 
energy distributions. Although there do appear to be separate core and halo 
components in model IV at t = 1 and t = 2, this distinction disappears as 
the bar develops. By t = 5 the energy distribution is nearly uniform over 
a substantial range. The similarity between the energy distributions at t = 4 
and t = 5 indicates that the latter distribution has essentially reached an 
equilibrium state. Such is not the case for the energy distribution of model V 
shown in figure 48. It is clear from the kinetic energy ratios that this exam- 
ple has not yet reached equilibrium at t = 3. Indeed the energy distribution 
at t = 2 is substantially different from the one at t = 3. Although the 
bar is well developed, several more rotations appear to be needed to achieve 
equilibrium. 

From figures 39 through 43 it is sometimes difficult to determine whether 
a bar is forming in the central region of the system. A convenient method for 
obtaining more quantitative information on nonaxi symmetric behavior is to plot 
the projected surface mass density for a number of cylindrical shells with dif- 
ferent radii as a function of longitude. This is done in figure 50 for models I 
to V. It can be seen that at t = 3 no barlike nonaxi symmetry is present for 
models I to III, whereas for models IV and V a pronounced bar structure is 
observed. 

N-body calculations and other data used by Ostriker and Peebles (ref. 8) 
indicate that galaxies with T cir/I w l > 0.14 are unstable to the bar-forming 
mode. The empirical evidence for this criterion has been gathered for systems 
near equilibrium. For the present calculations it makes more sense to inquire 
about the value of T c i r /|w] after the collapse phase rather than at t = 0, 
especially since the first evidence of bar formation seems to appear during the 
second collapse. The ratios T c i r /|w| for models I, II, IV, and v (figs. 46 
and 47) are either well below or well above the value 0.14 and they behave 
accordingly. The ratio for model III reaches a low of 0.13 at t = 0.75, and 
subsequently oscillates between 0.14 and 0.15. Nevertheless, this system 
remains axially symmetric as indicated most convincingly by figure 50. This 
apparent disagreement could easily be due to the crudeness of the empirical 
criterion or to slight numerical errors in the computation of the ratio. The 
most likely cause of such numerical errors is an underestimate of the potential 
energy. This occurs because the interaction potential used in these calcula- 
tions is cut off (and thus much lower than the actual Newtonian interaction 
potential) for distances less than one cell dimension (ref. 24) . This effect 
becomes more pronounced for systems with high central mass concentrations. 

Such an effect may also be responsible for the slightly high ratio obtained in 
the three-dimensional disk-core case. 
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Detailed Features 


We now turn to the more detailed diagnostics of these simulations. The 
volume density P is shown in figure 51 as a function of the spherical radius 
r for three values of the colatitude 0. This density represents an azimuthal 
average for a ring with radial thickness of 1 kpc and angular width of 10° in 
colatitude. The data points at r = 10 kpc for t = 0 appear to be a factor 
of 2 too small because the initial uniform density configuration fills only half 
of the ring centered at 10 kpc. 

For models I to III, the final density plots display the same separation 
into a core and a halo that was evident in the energy distributions. In the 
density plots, this property manifests itself as a change in slope between 
6 and 8 kpc. As before, this distinction becomes less sharp in the more rap- 
idly rotating example. In the case of model I, and to a lesser extent for 
model II, it is tempting to describe the density as the sum of two exponentials. 
This feature is displayed more clearly in plots of the density projected onto 
the equatorial plane as shown in figures 52 and 53. 

In the two bar models, the azimuthally averaged density in the equatorial 
plane maintains an exponential form up to at least 12 kpc. At higher inclina- 
tions, model V departs sooner from the exponential shape. This model, however, 
has not yet equilibrated. The exponential density profiles observed here have 
also habitually resulted from the evolution of flat galaxies into bars (ref. 4) . 

A different perspective of the densities is offered by figures 54 and 55, 
which show the variation in the vertical direction. Here, too, are the signs 
of a core-halo structure in the first three models and of a well-developed expo- 
nential profile in the last two. 

The rotation curves in the equatorial plane V(j> o = rw Q = \jrK r and the mean 

circular velocity are displayed in figures 56 and 57. The mean circular veloc- 
ity represents the mean velocity in a cylindrical shell covering all z. Note 
that heating, transfer of angular momentum to larger radii, and the increasing 
central density result in large reduction of the mean circular velocity when 
compared with the rotation curve ru) 0 . 

The mean rotation curves for various colatitudes are displayed in fig- 
ures 58 and 59. The azimuthal velocities shown there were averaged in the same 
manner as the densities given in figure 51. Model I, of course, has no net 
rotation. Models II and III have rotation curves which are linear for the 
inner 4 or 5 kpc and which turn over near 8 kpc. The location, shape, and 
magnitude of the peak in these two sets of rotation curves resemble closely 
those obtained by Gott (ref. 17). At both 0 = 90° and 0 = 45°, the present 
peak rotational velocities are within 10 km-sec -1 of those obtained by Gott. 

The reader who wishes to make his own comparison with figure 4 of Gott (ref. 17) 
should note that the virial velocity to be used there is 227 km-sec -1 and recall 
that the initial radius is 10 kpc. 

The final rotation curve for model IV is substantially different from the 
earlier axisymmetric work. The region of solid body rotation here extends for 
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only 1 or 2 kpc from the center rather than for a half dozen. Moreover, the 
peak in the rotational velocity occurs near 11 kpc rather than 6, and it is at 
least 20 km-sec -1 lower. None of this is at all surprising, considering the 
vastly different shape assumed by this model when the nonaxisymmetric forces 
are unleashed. On the other hand, the t = 2 version of the model IV rotation 
curve is quite similar to Gott’s result. After only 2 rotations the bar is 
still in its infancy while the collapse phase is nearly completed. As Gott 
noted, the violent relaxation of the axially symmetric collapse does tend to 
produce an extended region of solid body rotation. A comparison between the 
model IV rotation curves at t = 2 and t = 5 illustrates that the violent 
relaxation of the bar-forming stage produces quite a different result. The 
other bar system, model V, appears to be tending toward a rotation curve with 
properties similar to that of model IV. Given the additional integrals of 
motion that are likely to be present in a bar system, the theoretical interpre- 
tation of the last two rotation curves in terms of violent relaxation presents 
more difficulties. Note here that the rotation rates of the bars in models IV 
and V are 16 km-sec -1 -kpc -1 and 14 km-sec -1 -kpc -1 , respectively. 

Some of the data obtained for the velocity dispersions are presented in 
figures 60 to 65. The radial velocity dispersions for models I to III are again 
similar to those of Gott. They indicate a small isothermal region in the center 
and larger dispersions perpendicular to the equatorial plane than in it. Curi- 
ously, the present velocity dispersions have the same local maxima, for example, 
near r = 7 for model I, that were exhibited by Gott's results. 

These local maxima also appear in the vertically averaged velocity disper- 
sions, which are a cruder but more statistically reliable measure than the one 
presented here. For model IV, a comparison of the various velocity dispersions 
at t = 2 and at t = 5 gives some indication of the difference between the 
velocities produced by the collapse itself and those that result from the bar 
formation. At least in the azimuthally averaged sense displayed here, the bar 
is seen to have relatively little effect upon the velocity dispersions in the 
vertical direction while increasing those in the horizontal direction by 30 to 
40 percent near the center. 

We conclude this section with a discussion of some crucial observable 
quantities. Recent measurements of rotational velocities in ellipticals by 
Illingworth (ref. 21), Peterson (ref. 22), and others suggest a substantial 
disagreement between the relatively high rotational velocities that arise in 
the collisionless models which remain oblate, that is the models of Gott 
(refs. 17 and 18) and Wilson (ref. 35), and the surprisingly low rotational 
velocities which have actually been observed. Models I to III, which behave 
just like the corresponding models of Gott, share this deficiency. 

The bar models, on the other hand, offer some hope, as has been noted by 
Miller and Smith (ref. 20) and by Illingworth (ref. 21). We provide here some 
pertinent data on model IV, a case for which the bar is nearly equilibrated. 

The projected isodensity contours shown in figure 66 are especially illuminat- 
ing. There we have established a coordinate system relative to the rotating 
bar which retains the z-axis as the rotation axis but which aligns the y-axis 
with the largest axis of the bar (and, accordingly, places the x-axis along the 
shortest axis). The line-of-sight direction is thus specified by the usual 
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spherical angles. For instance, a viewing direction along the rotation axis 
has 0=0° while a line of sight in the equatorial plane is specified by 
0 = 90°. Moreover, a viewing direction which sees the bar broadside is speci- 
fied by 4> = 0° and one which sees the bar along the long axis is <J> = 90°. 

As noted earlier, if the long axis of the bar is taken to have unit length, then 
the vertical axis has length 0.54 while the short axis in the equatorial plane 
has length 0.39. 

The six views shown in figure 66 certainly give the impression of an ellip- 
tical of high eccentricity. The apparent ellipticity for these and a number of 
other viewing directions is listed in table II. The definition e = 1 - b/a, 


TABLE II.- ELLIPTICITY AND VELOCITY RATIO AS A FUNCTION OF 
VIEWING ANGLE FOR MODEL IV 


e, 

deg 

4>, 

deg 

e 

VmAMO) 
for slit along 
major axis 

Vm/o v (0) 
for slit in 
equatorial plane 

30 

0 


0.61 

0.61 


30 


.54 

.63 


45 


.47 

.64 


60 


.35 

.66 


90 

.71 

.10 

.65 

45 

0 

.68 

.92 

.92 


30 

.67 

.80 

.91 


45 

.66 

.67 

.89 


60 

.64 

.49 

.87 


90 

.65 

.10 

.83 

60 

0 

.65 

1 .20 

1.20 


30 

.62 

1 .07 

1.11 


45 

.58 

.89 

1 .07 


60 

.57 

.71 

1 .01 


90 

.54 

.10 

.97 

90 

0 

.61 

1 .48 

1 .48 


30 

.56 

1.29 

1 .29 


45 

.51 

1 .15 

1 .15 


60 

.28 

1 .07 

1 .07 


90 

.29 

.10 

1 .03 


where a and b are the measured long and short axes, is used. The values of 
e listed in table II are the average of the ellipticities for the two innermost 
isodensity contours. For viewing directions near the rotation axis, the ellip- 
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ticity is rather insensitive to the azimuthal orientation. When viewed from 
the equatorial plane, however, there is a substantial dependence on <j>. This 
effect is more dramatic for the thin bar which arose in model IV than it would 
be for a flat bar, one for which the vertical axis is the smallest. In the 
present case the projection effects also include a sudden switch of the apparent 
long axis. 

An important consequence of the triaxial nature of this bar is that the 
projected long axis of the bar can be quite misaligned with the equatorial plane. 
Figure 66 provides several vivid examples. All but the <J> = 90° and <j) = 0° 
cases have the apparent major axis substantially misaligned with the equatorial 
plane. This can have a sizable effect upon the observed velocities if the mea- 
surements are made only with the slit aligned along the apparent major axis. We 
have imposed a rectangular slit upon many projected views of model IV. This 
slit was 6 kpc wide and 32 kpc long and was centered upon the center of the sys- 
tem. The average line-of-sight velocities were computed in 1 kpc x 6 kpc por- 
tions along the slit. The peak rotational velocity V m and the central line- 
of-sight velocity dispersion a v (0) were then computed. These measurements 
were performed for slits oriented both along the apparent major axis and perpen- 
dicular to the rotation axis. The results for the ratio U = V m /a v (0) are 
given in table II. The slit oriented perpendicular to the rotation axis picks 
up the maximum rotation. As the viewing direction varies, the changes in the 
average line-of-sight velocities reflect not only the geometric effects which 
occur even in oblate models, but also the noncircular streaming motions along 
the bar. When just these two effects are considered, the ratio U still seems 
to be fairly high for most orientations. But when the shifted orientation of 
the apparent major axis is also taken into account, the ratio u seems much 
closer to observed values. Note that the ellipticities for these orientations 
of model IV are all on the high side so that a U « 0.50 to 0.70 appears con- 
sistent with reported measurements. Nevertheless, a rough estimate suggests 
that a randomly chosen viewing direction will have U < e about 40 percent of 
the time. 

The same six projections are shown for model V in figure 67. Weak spiral 
features are still evident at t = 3, indicating that the bar has not yet 
reached equilibrium. Note again how thin the bar is. 

Figure 68 displays polar and equatorial projections for models I to III. 
These show both the increased flattening and the decreased central concentration 
caused by increasing amounts of initial rotation. 


DISCUSSION AND CONCLUSIONS 

Comparison of the evolution of an infinitesimally thin stellar disk with 
a three-dimensional stellar disk for an initially exponential radial density 
variation shows that the primary effect of finite thickness is a reduction in 
the initial growth rate of the spiral or bar instability. Vandervoort's 
(ref. 36) formula for the thickness correction to the dispersion relation for 
an infinitesimally thin disk suggests that the effective Q for the three- 
dimensional system is the range 1.3 to 1.4 compared with roughly 1.0 for the 
infinitesimally thin disk. Nevertheless, after a time of three rotational 
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periods the systems have evolved to essentially the same general structure. 

Thus, two-dimensional disk models of galaxies are appropriate for the investi- 
gation of various global, spiral structure, or instability problems associated 
with flat galaxies. Adding a fully self-consistent central, nonrotating, expo- 
nential core/halo component has a considerable stabilizing effect on the system. 
The effective Q for this second three-dimensional system is in the range 1.5 
to 1 .6 from r = 5 kpc outward and considerably higher inward. While for the 
two-dimensional system an outer spiral and central bar structure quickly forms, 
the three-dimensional system's most apparent instability is a weak spiral struc- 
ture evolving from Jeans' type instabilities in the outer regions of the disk. 
The spiral structure for the three-dimensional disk-core/halo system is tending 
to dissipate by the time bar formation is well in progress, as determined from 
the azimuthal density variation. Thus, even for systems with large core/halo 
components, two-dimensional disk models appear sufficient to simulate evolution 
of spiral and other global structures of disklike galaxies. The primary effect 
of performing the simulation in three dimensions is a reduction in the growth 
rate of global instabilities. 

The range of initial conditions for this study of the fully three- 
dimensional, collisionless collapse of uniform spheres of stars covered some 
models which remained axisymmetric and others which developed bar-shaped struc- 
tures. The first four models had been treated previously by Gott in a manner 
which precluded the development of nonaxisymmetric features. Our results for 
models 1 to III showed that where large-scale axial symmetry was retained, even 
the detailed features were well represented by the strictly axisymmetric calcu- 
lations. The energy distributions, density profiles, systematic velocities, 
and random velocities were all similar to those obtained by Gott. Moreover, 
no significant redistribution of angular momentum occurred, even though per- 
mitted in principle by the numerical treatment. Nevertheless, the present 
results for models IV and V illustrate the limitations of the axisymmetric 
calculation. These limits had been predicted reliably by Gott himself using 
the Ostriker-Peebles criterion. For model III, the ratio of kinetic energy in 
rotation to total potential energy was too close to the Ostriker-Peebles sta- 
bility limit of 0.14 to make a definitive judgment. The present results sug- 
gest once again that the natural outcome of strong nonaxisymmetric instabilities 
is bar formation. 

The flattest oblate system obtained here was model III, which resembles 
the appearance of an E2 galaxy. The bar instability of model IV prevented the 
development of an oblate system as flat as an E5, such as Gott had done with 
his constrained model. On the other hand, bar formation led to even flatter, 
albeit triaxial, systems. Model V, which can appear as flat as an E8, is the 
more extreme example. 

The bar-forming models display a distinct separation between the collapse 
and the bar phases. The bars appear to form during or after the second collapse 
along the rotation axis. Prior to this time the evolution is dominated by axi- 
ally symmetric collapse dynamics. Model IV provides an excellent example of 
this since a good deal of similarity prevails between the properties of the 
present model IV at t = 2 and the final form of the axially constrained ver- 
sion. This same comparison made at t = 5 suggests the following effects of 
the bar upon the detailed properties: (1) the distinct core and halo components 
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produced by the collapse are obliterated by the bar; (2) the density is more 
clearly represented as a combination of two exponentials; (3) the central region 
of solid body rotation is substantially reduced; (4) the dispersions of the hori- 
zontal velocity components are increased by 30 to 40 percent while the disper- 
sions of the vertical component are nearly unaffected. 

Given the recent observational data on the apparent low systematic veloci- 
ties in ellipticals, the bar models seem much more viable for galaxy formation 
than the oblate ones. We have explored extensively the equilibrium properties 
of only one bar model. Although there are many viewing directions which are 
reconcilable with the observations (at least observations made solely along the 
apparent major axis) , there are uncomfortably many others which would suggest 
high rotation. No doubt, a more satisfactory model could be constructed by 
using different initial conditions than those employed in model IV. Certainly 
less rotation should be applied and perhaps even slightly greater initial veloc- 
ity dispersion. 

It is clear that there is one observational test which even a finely tuned 
triaxial model must pass. The data in table II indicate how sensitive the mea- 
sured rotation of the bar models can be to the slit orientation. Thus, for many 
real triaxial systems of the type modeled by the final two examples, there ought 
to be some slit orientations which produce substantially greater measured rota- 
tion than a slit orientation corresponding to the apparent major axis. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
March 6, 1979 
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APPENDIX 


COMPUTER PROGRAM FOR GENERATING THE THREE-DIMENSIONAL GRAVITATIONAL 
POTENTIAL DISTRIBUTION OF ISOLATED GALAXIES 
Mathematical Summary 

The scaled gravitational potential at the center of cell (x,y,z) is defined 
by the triple summation over the three-dimensional array of cells 


2n-l 2n-l 2h-1 

^ ^i, j , k H i-x, j-y ,k-z ( A1 ) 

i=0 j=0 k=0 

where 

H i,j,k = (i 2 + + k 2 ) ^ 2 for i + j + k f 0 

H 0 , 0,0 = 1 


***>*■ = A Z 


and Pi f j f k is the mass density in cell (i,j,k). Because direct summation 
is much too time consuming to be practical, the triple summation is evaluated 
by the convolution method using fast Fourier transforms (ref. 2) . That is, the 
Fourier transform of the potential equals the product of the Fourier transforms 
of P and H: 


*C,n,C ■ p C,n,^,n,c 


(A2) 


The gravitational potential 4> x ,y,z is obtained by taking the inverse Fourier 
transform of equation (A2) . Rather than the usual complex Fourier series, a 
real expansion is used here. For example, the Fourier transform of the density 
Px,y,z is given by 


2h-l 2n-l 2n-l 


p £,n,C = X X X C(X,n> c(y,n) c(z ' h > p x,y,z f(Ax,n) 


z=0 y=0 x=0 


x f(n,y,n) f(C,z,h) 


(A3) 
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where 


f (£,x,n) 


cos 

,sin 


(Tr^x/n) 
tt - n) x/n 


0 S 5 S n 
n < ^ 2n - 1 


c(x,n) 



if x = 0 or x = n 
Otherwise 


The symbols n and h define the n x n x h active array and also the 
(2n) x (2n) x (2h) larger array over which the Fourier transform must be taken 
so that the potential for an isolated galaxy is obtained (see fig. 69). Note 
that the density may be nonzero only in the smaller n x n x h array. Because 
of the symmetry of H X/ y /Z , the Fourier transform can be obtained by a 

finite cosine transform 


h n n 

= X Z- zE c2(X ' n) c2(y,n) c2 < z ' h > H x,y ,z 
2=0 y=0 x=0 


x cos (7r^x/n) cos (rrriy/n) cos (Tr^z/h) 


(A4) 


for 0 i CfH = n and 0 S x, % h, and 


H C+n,n,C H C+n,n+n,C H £+n,n f C+h ~ H £+n,n+n,C+h 
= H C,n+n,C = ^,n+n,C+h = ^, n , 5 +h = 


The next step in obtaining the potential is to multiply P£,n,£ by H £,ri,£ 
to obtain < t > £,n,£ ( e< 3« (A2)). 

The gravitational potential for an isolated galaxy correctly defined over 
the n x n x h array is obtained by the Fourier synthesis 


$x,y,z 


1 

N 3 


2h-l 

I 


C=0 


2n-l 


I 

"I =f 


2n-l 

C=o 


f(H,y,n) f (C ,z,h) 


(A5) 
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Note also, that since 


H C,u,c - H S,s,n = = 


different permutations of the same set of indices need not be stored. Thus, 
the transformed Green's function can be converted to a one-dimensional array 


H S,n,c 


= F 


n 


where different permutations of are stored in the same location n 

given by 


n 



n 

i) + - (n - i ) + e 
2 


= C(5 - 1) (2£ - 1)/12 + - l)/4 + n (n - l)/2 + c 


Computer Program Subroutine Which Uses Only Core Storage 

At the end of this appendix, FORTRAN listing I is given of computer program 
GETPHI which may be used to obtain the potential <t> by use of a (2n) x (2n) x h 
array of cells. The variable I2A defines the x and y dimensions, and I3A 
the z dimension of the array used for the potential calculations. When the 
subroutine GETPHI is called, RHO(I,J,K) contains the mass density and GETPHI 
places the values of the corresponding gravitational potential in RHO(I,J,K). 

The subroutine FTRANS (I, I2B) , written by R. Hockney (ref. 37), performs a finite 
Fourier analysis or synthesis on the common input array Z and places the 
result in the common output array Y. The subroutine performs a cosine analysis 
for I = 2, a periodic analysis for 1=3, and a periodic synthesis for 1=4. 
The subroutine GETSET (I, I2B) initializes FTRANS and is called every time the 
arguments of FTRANS (I , I2B) are changed. The Fourier transform is cal- 

culated on an (n + 1 ) x (n + 1 ) x (h + 1 ) array only the first time that the 
subroutine is called and is kept in storage for subsequent use. 

The Fourier transform of P x ,y,z * n ^e x-direction is generated by obtain- 
ing the partial transform P%,y r z f° r 0 S £ ^ 2n - 1 , 0 S y £ n - 1, and 

0 S z £ h - 1 . The transform P^,y, z is zero outside of this region because 
Px,y,z i s nonzero only over the n x n x h active array. Next, the Fourier 
transform of P£ f y, z is performed in the y-direction obtaining the x-y partial 
transform P£ f q, z r for 0£££2n-l, 0 = U = 2n - 1 , and 0 S z ^ h - 1 . 

Since ,n , z i s zero for h ^ z S 2h - 1 , by use of one-dimensional arrays Y 
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and Z the Fourier transform of can be taken in the z-direction to 

obtain the total transform P£ n ,T,~ f° r 0 = ? = 2h - 1 . Next, P£ f n,£ is 
multiplied by H£ 7 n,£ to obtain and 4:116 i nverse Fourier transform 

is performed in the z-direction. The resulting partial x-y transform ^ z 
is placed in the 2n x 2n k h RHO(I,J,K) array for 0 = K = 2n - 1 , ’ 

0 S n = 2n - 1 , and 0 = z = h - 1 with values for h = z = 2h - 1 discarded. 
(The use of these one -dimensional arrays was first presented in ref. 23 for a 
two-dimensional potential solver.) Next, the inverse Fourier transform of 
is generated in the y-direction to obtain the x partial transform 

$£,y, z for 0 = £ = 2n - 1 , 0=y=n-l, and 0 = z = h - 1. The final 
step is to perform the inverse Fourier transform in the x-direction for 
0 = y i n - 1 and 0 = z = h - 1 to yield the correct gravitational potential 
< t > x,y,z for an isolated galaxy over the n x n. x h array., 

Overlayed Computer Program Which Uses Core and Disk Storage 

Use of program GETPHI, with the 64 x 64 x 16 active density /potential array 
used in this report would have necessitated the dimensioning of the RHO array 
at 128 x 128 x 16 and the H array at 65 x 65 x 17. Since such large dimen- 
sions would have precluded use of the CDC 6600 computer, program GETPHI was 
modified to include use of overlayed programs and disk storage resulting in a 
maximum core storage requirement at any one time equal to about five-fourths 
of the active array. The listing of this program (listing II) includes (1) a 
section of an initializing overlay in which relevant constants are computed, 

(2) a section of the star advancing overlay in which "chunks" of the density 
array are written on appropriate disk files, (3) another section of the star 
advancing overlay in which chunks of the computed potential array are read 
from disk files, (4) the GETH overlay which computes H, and (5) the GETPHI 
overlay which computes the potential array from the density array. Figures 70 
through 74 are presented to facilitate description of the overlayed program. 

For clarity, figures 69 through 74 are drawn for an active array dimensioned 
nxnxh=8x8x4; table III compares the array dimensions of these figures 
with those of the overlayed program. 

The method used is the alignment in the direction of transformation of four 
identical arrays named RH01 , RH02, RH03, and RH04, each of which is dimensioned 
(n/2) x (n/2) x h within the GETPHI overlay. (See figs. 70 and 73.) The active 
array is dimensioned as the PHI array within the initializing and star advancing 
overlays (see figs. 69 and 71) but is not dimensioned within the GETPHI overlay. 
As figure 70 suggests, the chunks RH01 , RH02, RH03, and RH04 may be visualized 
as forming either a row or a column of the lower half (0 = z = h - 1 ) of the 
extended array. Switching the lineup to a different row or column is accomp- 
lished by storing the array associated with each chunk location on a separate 
file; these eight files are also indicated in figure 70. 

As shown in figure 71 , one chunk size array named 01 is dimensioned in 
the initializing and star advancing overlays. Chunks of the active array are 
transferred between the PHI array of these overlays and the arrays RH01 , RH02, 
RH03, and RH04 of the GETPHI overlay via "do loop" transfer to/from the 01 
array and storage on files 1, 2, 5, and 6. 
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TABLE III.- ARRAY DIMENSIONS FOR PROGRAM OF LISTING II 
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At the beginning of a program run, the GETH overlay computes H in the 
(n + 1 ) x (n + 1 ) x (h + 1 ) H array in the same manner as listing I. All of 

except for two boundary planes of elements (£ = n, 0 = h = n, 0 = C = h, 
and 'o ^ £, = n, n = n, 0 = £ = h) is then transferred in portions via "do loop" 
to the (n/2) x (n/2) x (h + 1) HH array from which it is^written on disk 
file 9 (see fig. 72). Elements of one boundary plane of H£ fT1f £ (K = n, 

0 = n=n,0S£Sh) are transferred to the (n + 1 ) x (h + 1) HN21 array 
which is in common with the GETPHI overlay; the C-h transpose of that boundary 
plane is equal to the other boundary plane (0=£=n,n=n, 0 = ? = h) because 
of the symmetry of H across the £=n diagonal plane. During each poten- 
tial solution the portions of H on file 9 are read sequentially into an 
(n/2) x (n/2) x (h + 1) HH array of the GETPHI overlay from which H elements, 
along with those in the HN21 array, are multiplied with p. This sequence 
(listed in table IV) utilizes the symmetry and periodicity of H (eq. (A4)) to 
provide a full set of (2n) x (2n) x (2h) H elements to the GETPHI overlay 
in a manner which minimizes the reading of file 9. 

TABLE IV.- STORAGE OF THE FOURIER TRANSFORMED GREEN'S FUNCTION H 

ON DISK FILE 9 
[program of listing if] 


Record no. 
of file 9 

Storage sequence within 
GETH overlay 

(a) 

Use sequence within GETPHI 
overlay 
(b) 

1 

A 

(1,1), 0,3) 

2 

B 

(1,2), (1,4), (3,2), (3,4) 

3 

A 

(3,1), (3,3) 

4 

C 

(2,1), (2,3) 

5 

D 

(2,2), (2,4), (4,2), (4,4) 

6 

C 

(4,1), (4,3) 


a Within the GETH overlay, this is the location in the H array (as 
designated by letters A-D of fig. 72) from which "do loop" transfer is 
made to the HH array followed by writing on the indicated record of 
disk file 9. 

^Following reading of the indicated record of disk file 9 into the 
HH array within the GETPHI overlay, this is the sequence of locations 
in the extended PHI array (as designated by "chunks" (IROW, JCOLUMN) of 
fig. 70) upon which z-direction one-dimensional array operations are 
performed. These operations include multiplications by H, the appro- 
priate portion of which is now contained in the HH array. This method 
minimizes reading of file 9 by using the periodicity and symmetry of H. 
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The GETPHI overlay consists of subroutines ANLX ( JCOLUMN) , ANLSYN (IROW) , 
and SYNX (JCOLUMN) which dimension in common the arrays HH, HN21 , RH01 , RH02, 
RH03, and RH04, as pictured in figure 73. Figure 74 indicates the lineup of 
chunks associated with each call to a subroutine. The potential solution is 
mathematically identical with that described for listing I. Calling ANLX(l) 
and ANLX (2) performs the Fourier transform of Px,y,z x-direction to 

form P£, y , z . Calling ANLSYN (1 ) , ANLSYN (3) , ANLSYN (2) , and ANLSYN (4) in 
sequence performs the following: (1) a Fourier transform of P£,y, z in the 

y- and z-directions to form P£,n,i;» ( 2 ) multiplication with to form 

and (3) the inverse Fourier transform of 4>£,n,C t * le z ~ an< 3 
y-directions to form $£,y,z* Calling SYNX(l) and SYNX(2) performs the 
inverse Fourier transform of $^,y,z i n x-direction to form 4>x,y,z* 

The GETPHI overlay is outlined in more detail in table V. 


Comparison of the Two Computer Programs 

The overlayed program (listing II) is preferable to the original program 
(listing I) because the addition of some peripheral processing time and a 
small increase in central processing time are much more than compensated by a 
75-percent decrease in the required core storage. The maximum number of active 
array elements dimensionable on the CDC 6600 with the original program is 16384 
(e.g., 32 * 32 x 16) and with the overlayed program, 65536 (e.g., 64 x 64 x 16); 
the latter program can have other potentially useful active array dimensions 
of 32 x 32 x 8, 32 x 32 x 16, and 32 x 32 x 32. Solution of the 64 x 64 x 16 
active array by the CDC 6600 requires about 300 000 (octal) words of core stor- 
age and with H already computed takes about 75 seconds of central processing 
time. 
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TABLE V.- OUTLINE OF THE GETPHI OVERLAY OF THE PROGRAM OF LISTING II 

Refer to figure 74 for orientation of arrays RHOl , RH02, RH03, 
and RH04 and to figure 70 for file numbers corresponding to 
the locations of these arrays 


Listing II 
line nos. 


A. CALL ANLX(l) (fig. 74(a)) 

1. Read files 1 and 2 into RHOl and RH02, respectively 299-302 

2. Set RH03 = RH04 = 0 316-317 

3. Perform Fourier transform in x-direction over RHOl, RH02, 

RH03, and RH04: P x ,y,z ^ P£,y,z 310-323 

4. Write RHOl, RH02, RH03, and RH04 onto files 1, 2, 3, 

and 4, respectively 325-332 

B. CALL ANLX ( 2) (fig. 74(b)) 

1 . Read files 5 and 6 into RHOl and RH02, respectively 305-308 

2. Same as steps A. 2 and A. 3 

3. Write RHOl, RH02, RH03, and RH04 onto files 5, 6, 7, 

and 8, respectively 335-342 

C. CALL ANLSYN(1 ) (fig. 74(c)) 

1. Read files 1 and 5 into RHOl and RH02, respectively ..... 353-356 

2. Set RH03 = RH04 = 0 382-383 

3. Perform Fourier transform in y-direction over RHOl, RH02, 

RH03, and RH04: P£, y , z -► P£,n,z 376-389 

4. Read record 1 of file 9 into HH 393 

5. For each one-dimensional array in z-direction of which 

RHOl is composed 394-407, 


477-483 

a. Transfer to one-dimensional array Z, dimensioned at 

least 2h + 1 

b. Set Z = 0 for z s? h 

c. Perform Fourier transform in z-direction over Z for 

0 S z = 2h - 1 with the result appearing in one- 
dimensional array Y: P£,p,z P£,r|,£ 

d. Multiply Y by to form = P^, n , c %, n , ? 

e. Perform inverse Fourier transform in z-direction over Y 


and store^result for 0 ^ z £ h - 1 in RHOl : 
h, n,z 

6. Repeat step C.5 for RH03 426-454, 

471-483 

7. Read record 2 of file 9 into HH 410 

8. Repeat step C.5 for RH02 and RH04 411-424, 

456-469, 

477-483 

9. Perform inverse Fourier transform in y-direction over 

RHOl, RH02, RH03, and RH04: $^,n,z ^ $£,y, z 486-497 


10. Write RHOl and RH02 onto files 1 and 5, respectively .... 500-503 
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TABLE V.- Concluded 


Listing II 
line nos. 

D. CALL ANLSYN (3) (fig. 74(e)) 

1 . Read files 3 and 7 into RHOl and RH02,. respectively 365-368 

2. Same as steps C.2 to C.9 except for sequencing of reading 

tape 9 into HH and the z-directional operations. 

Table IV details this sequencing. 

3. Write RH01 and RH02 onto files 3 and 7, respectively .... 512-515 

E. CALL ANLSYN (2) (fig. 74(d)) 

1 . Same as step D except that files 2 and 6 correspond to 
RH01 and RH02, respectively, for read and write 
operations 

F. CALL ANLSYN (4) (fig. 74(f)) 

1 . Same as step D except that files 4 and 8 correspond to 
RH01 and RH02, respectively, for read and write 
operations 


G. CALL SYNX(l) (fig. 74(a)) 

1. Read files 1, 2, 3, and 4 into RHOl , RH02, RH03, 

and RH04 , respectively 530-537 

2. Perform inverse Fourier transform in x-direction over 

RHOl, RH02, RH03, and RH04: <t>£, y , z ^ 4> x ,y,z 550-560 


3. Write RHOl and RH02 onto files 1 and 2, respectively .... 562-565 

H. CALL SYNX (2) (fig. 74(b)) 

1. Read files 5, 6, 7, and 8 into RHOl, RH02, RH03, 

and RH04 , respectively 540-547 

2. Same as step G.2 

3. Write RHOl and RH02 onto files 5 and 6, respectively .... 568-571 
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LISTING I.- SUBROUTINE GETPHI FOR CALCULATING THE THREE-DIMENSIONAL 
GRAVITATIONAL POTENTIAL USING ONLY CORE STORAGE 


subroutine getphi 

COMMON Z( 1 02 5 ) . Y( 1025 ) .RH0I64 .64 . I 6 1 • 1 2A . 1 3A . I TEST 
DIMENSION HI33.33.17) 

IF ( I TEST.EQ.O ) GO TO 11 

1 TEST=0 

128= 12A-I 

N = 2** 1 2 A 

N02=n/2 

N2l = N02 + 1 

I 38= I 3A-1 

NH=2** 1 3A 

NH02=NH/2 

NH21 =NH02+1 

RNI = 1 ./ <N*N*NH ) 

DO I K = 1 .NH21 
DO 1 J=1,N21 

DO 1 1=1, N21 

R I = ( K- 1 ) * ( K — l )+( J-l >*(J-1 )•*•< 1-1 )*( 1-1 ) 

IFIRI .LT.l . ) R I = 1 • 

H ( I , J.K )=RNI /SORT (R 1 ) 

1 CONTINUE 

CAUL GETSET (2 . 12B ) 

DO 2 <=1 »NH21 
DO 2 J= 1 .N21 
DO 3 1=1. N2 1 

3 21 I >=H< 1 , J.K ) 

CALL FTRANS ( 2 . 128 ) 

DO 4 1=1. N2 1 

4 HI | , J.K ) = Y ( I ) 

2 CONTINUE 

DO 5 K = 1 ,NH21 
00 5 1=1. N2 1 
OO 6 J= 1 .N21 

6 Z< J)=H( I .J.K ) 

CALL FTRANSI2. I 26 ) 

00 7 J=1 .N21 

7 HI I . J.K >=Y< J 1 

5 CONTINUE 

CALL GETSET 12 . 1 38 ) 

DO 10 J=1 ,N21 
DO 10 1=1 .N21 
DO 8 K = 1 , NH2 1 

8 2 IK ) =H I I .J.K) 

CALL FTRANS 12. 138) 

DO 9 K = 1 , NH2 1 

9 HI l .J.K ) = Y I K ) 

10 CONTINUE 

11 CONTINUE 
WRITEI6.43) 

43 FQRMATI10H HII.J.K)) 

DO 42 K = 1 .NH21 
OO 42 J=1 .N21 
WR I TE 16.41 ) J.K 

WRITE16.40) (HI I . J.K ). 1 = 1 ,N21 ) 

41 FORMAT I 1 4H 1 = 1. N21 J=13.5H K=I3) 

40 FORMAT 1 2H 8E16.8) 

42 CONTINUE 

CALL GETSET (3, I2A ) 

OO 14 K = 1 ,NH02 
00 14 J=1 .N02 
DO 12 I =1 .N 

12 Z( 1 ) =RHO I I .J.K ) 

CALL FTRANS I 3 . I2A ) 

DQ 13 I =1 ,N 

13 RHO I I ,J,K)=Y( I ) 

14 CONTINUE 
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DO 17 KM « NH02 
DO 17 1 =1 *N 
DO 15 J=1 *N 

15 Z( J)=RHO( I .J.K ) 

CALL FTRANS (3. 1 2 A ) 

DO 16 J=1 .N 

16 RHO< I .J.K ) = V ( J ) 

17 CONTINUE 

DO 20 1 = 1 nN 

DO 20 J= 1 . N 
DO 18 K=1 .NH02 
Z IK ) =RHO ( I .J.K) 
ie Z (K4-NH02 ) =0, 

CALL GETSET (3* 13A ) 

CALL FTRANS (3. 13A ) 

IF ( ! .GT .N21 .AND*J*LE«N2l ) GO TO 22 
IF ( I .LE.N21 .AND#J.GT.N21 ) GO TO 24 
IF ( ! .GT *N21 .AND* J.GT.N21 ) GO TO 26 
DO 19 K = 1 « NH02 
Z(K)=Y<K>*H< I .J.K > 

19 Z (K+NH02 > = Y (K+NH02 >*H< I • J.K ) 

Z ( 1 ) =Y { 1 )*H( I • J . 1 ) 

Z ( NH2 1 ) = Y ( NH2 1 ) *H ( I • J • NH2 1 ) 

GO TO 21 

22 DO 23 K =2 • NH02 
Z(K)=Y(K)*H( I — N02 .J.K) 

23 Z IK+NH02 > = Y (K + NH02 )*H ( I — N02 • J *K ) 

Z( 1 ) =Y< 1 >*H< I -N02 * J • 1 ) 

Z ( NH2 1 > = Y < NH2 1 ) *H I I -N02 • J • NH2 1 > 

GO TO 21 

24 DO 25 K=2*NH02 

Z (K) =Y (K >*H( I ♦ J-N02.K > 

25 Z < K+NH02 )=Y(K>NH02 )*H< I . J-N02.K ) 

Z< 1 >=Y{ 1 )*H( I * J— N02 * 1 ) 

Z < NH2 1 )=Y (NH21 >*H ( 1 • J-N02.NH21 ) 

GO TO 21 

26 DO 27 K=2.Nh02 
ZIK)=Y<K)*H( 1-N02* J-N02.K) 

27 Z I K+NH02 > =Y (K+NH02 )*H< I -NQ2 • J-N02 • K 1 
Z < 1 > =Y { 1 ) *H ( I -N02 • J-NC2 . 1 ) 

Z (NH21 ) =Y (NH21 >*H ( I -NC2 * J-N02 .NH2 1 ) 
21 CONTINUE 

CALL GETSET (4. !3A ) 

CALL FTRANS (4 , I3A ) 

DO 28 KM *NH02 

28 RHO( I • J*K)=YCK ) 

20 CONTINUE 

CALL GETSET (4 , I2A ) 

DO 29 KM .NH02 
DO 29 J= 1 . N 
DO % 30 I M • N 

30 Z ( I ) =RHO < I « J .K ) 

CALL FTRANS (4 . I2A ) 

DO 31 I M .N 

31 RHO ( I . J.K )*Y ( I ) 

29 CONTINUE 

DO 32 KM . Nh02 
DO 32 1=1 .N02 
DO 33 J= 1 *N 

33 Z ( J)=RHO( I . J.K ) 

CALL FTRANS ( 4 . J2A ) 

DO 34 JM • N02 

34 RHO ( I . J.K )=Y< J) 

32 CONTINUE 
RETURN 
END 
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LISTING II.- OVERLAYS FOR CALCULATING THE THREE-D IMENS ION AL GRAVITATIONAL 

POTENTIAL USING CORE AND DISK STORAGE 

c the following is the section of an initializing overlay in which constants 001 


C RELATtO TO THE DIMENSIONS OF THE PHI ( DENS I T Y/POTEN T I A L ) ARRAY ARE COM- 002 

C PUTED- IT IS CALLED ONCE AT THE BEGINNING OF A PROGRAM RUN. IN THIS 003 

C LISTING THE VALUES OF 12 A* I 3 A AND THE DIMENSION AND LABELED COMMON 004 

C STATEMENTS ARE SET FOR AN ACTIVE PHI ARRAY DIMENSIONED 64 BY 64 BY 16* 005 

I 2 A = 7 006 

I 3 A = 5 007 

126= I 2 A - 1 008 

I 3B= I 3 A - 1 009 

N=2**I2A 010 

N02=N/2 01 1 

N2 1 =N02 + 1 012 

N04=N/4 013 

N34=N02+NO4 014 

NH=2**!3A 015 

NH02=NH/2 016 

NH2 1 =NH02+ 1 017 

C 018 

019 

020 


0 ************************************************************************** 021 
0 *************************************************************************** 022 
C THE FOLLOWING IS THE SECTION OF THE STAR ADVANCING OVERLAY IN WHICH CHUNKS 023 
C OF THE PHI ARRAY (CONTAINING THE DENSITY MESH) ARE WRITTEN ONTO DISK FILES 024 


C 1.2.5 AND 6 • THh STAR ADVANCING OVERLAY IS CALLED ONCL PER T I l>i(_ STEP. 025 

DIMENSION PM I (64 ,64 , 1 6 ) .Ol ( 32. J2 . I 6 > 026 

DO 520 K = 1 . NH02 027 

DO 520 J= 1 , N04 028 

00 520 I = 1 .NOA 0 29 

520 O I < I . J.K >=PH I < ! . J ,K ) 0 30 

WRITE! 1 ) 01 031 

REWIND 1 032 

DO 525 K = 1 . NH02 033 

DO 525 J= 1 . NOA 03A 

DO 525 1=1 .NOA 035 

525 01 < 1 , J.K )=PH I ( I .NOA+J.K ) 036 

WRITEC5) 01 037 

REWINO 5 038 

DO 530 K = 1 « NH02 039 

DO 530 J=1 .NOA OAO 

DO 530 1 =1 .NOA 0A1 

530 01 ( I . J.K >=PHI (NOA+ I . J.K I 0A2 

WRITEI2) 01 0A3 

REWIND 2 OAA 

DO 535 K=1.NH02 0A5 

DO 535 J= 1 .NOA 0A6 

DO 535 1=1 .NOA 047 

535 O I ( I . J • K ) =PH I < NOA+ I . NOA+ J . K ) 048 

WRITEI6) 01 049 

REWIND 6 050 

C 051 

C 052 

C 053 


q *****#******■*♦*****»***»**»»*#**-»■»**-»*# »********#*■*■■»******«**#♦»«*******»»* 05A 
C**********»»»-»**********»********-»* *********************** A**************** 055 
C THE FOLLOWING IS THE SlCTION OK THE STAR ADVANCING OVERLAY IN WHICH CHUNKS 056 
C OF THE PHI ARRAY (CONTAINING THE POTENTIAL MESH) ARE READ FROM DISK FILES 057 
C 1.2.5 AND 6. 058 
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DIMENSION PH I (64 . 64 . 16 ) .01 ( 32.32 . 1 6 ) 059 

READ ( 1 ) 01 060 

REWIND 1 061 

DO 30 K = 1 ,NH02 062 

DO 30 J=1,N04 063 

DO 30 1=1. N04 064 

30 PHI ( I . J.K )=OI C 1 . J.K ) 065 

REA0(5) 01 066 

REWIND 5 067 

DO 40 K = 1 ,NH02 068 

00 40 J= 1 . N04 069 

DO 40 |=1 . N04 070 

40 PH I ( I .N04 + J.K )=0! ( I . J.K ) 071 

READ ( 2 ) 01 072 

REWIND 2 073 

00 50 K= 1 . NH02 074 

DO 50 J= 1 . N04 075 

DO 50 1=1 . N04 076 

50 PHI (N04+ I , J.K 1 =01 ( ! . J.K ) 077 

READ (6) 01 078 

REWIND 6 079 

DO 60 K= 1 .NH02 080 

DO 60 J=1 . N04 081 

DO 60 1=1 ,N04 082 

60 PHI (N04+ I .N04 + J.K 1=0 I ( I . J.K 1 083 

084 

085 

C 086 

(;****■»******»*******************»#■»********#***♦*********#***«*♦*«*******♦** 087 

C****** **************************************************************** *»•»•**** 088 

C THE FOLLOWING IS THE GETH OVERLAY. WHICH COMPUTES AND STORES THE TRANS- 089 

C FORMED GREENS FUNCTION. IT IS CALLED ONCE AT THE 8EG INNING OF A PROGRAM 090 

C RUN. 091 

OVERLAY! 1FILE. 4.0 1 092 

PROGRAM GETH 093 

C THIS OVERLAY PERFORMS A COSINE ANALYSIS OF THE THREE-DIMENSIONAL GREENS 094 
C FUNCTION ARRAY. IT THEN WHITES CHUNKS OF THIS ARRAY ON DISK FILE 9 IN THE 095 
C ORDER IN WHICH THEY WILL HE READ INTO THt HH ARRAY DURING THE GETPHI 096 

C OVERLAY. VALUES FOR l=N/2+l AND J-N/2+ 1 ARE TRANSFERRED TO THE HN21 ARRAY 097 
C WHICH IS IN COMMON WITH THE GETPHI OVERLAY. 098 

COMMON/ALLCOM/N . NOE «N2 1 . N04 . N34 . NH. NH02 . NH21 .I2A.I2B.I3A.I3B 099 

COMMON/HN2 1 C0M/HN2 1 ( 65 .171 1 00 

COMMON 2(1025), Y ( 1 025 1 101 

DIMENSION H (65.65 . 1 7 1 .HH ( 32 « 32. 1 7 ) 102 

RN I = 1 . / (N*N*NH ) 103 

DO 1 K= 1 ,NH2 1 I 04 

DO 1 J= 1 ,N21 105 

DO 1 I = 1 . N21 106 

RI=(K-1)*(K-1)+(J-I)*(J-1)+(I-1)*(I-1) 107 

IF < R I .LT. 1 . 1 RI = 1 . 108 

H ( I , J.K 1 =RN I /SORT (R I 1 109 

1 CONTINUE 1 1 0 

CALL GETSET (2. I 28 1 111 

00 2 K=1 .NH21 1 12 

DO 2 J=1 .N21 113 

DO 3 I = 1 . N2 1 114 

3 Z( I )=H( 1 , J.K 1 115 
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CALL FTRANS (2 . I2B ) 116 

DO 4 l=l.N21 117 

4 H( I , J.K )=Y< I ) 118 

2 CONTINUE 1 19 

DO 5 K=1 .NH21 120 

DO 5 1=1 .N21 121 

DO 6 J=DN2I 122 

6 Z<J)=H( I .J.K) 123 

CALL FTRAnS<2. 126) 124 

DO 7 J= 1 .N21 125 

7 H< I . J.K )=Y< J ) 126 

5 CONTINUE 127 

CALL GETSETI2. I3B ) 128 

DO 10 J=1 ,N21 129 

DO 10 1=1 ,N21 1 30 

DO 8 K=1 • NH2 1 131 

a Z IKI=HI| , J.K > 132 

CALL FTRANS (2.1 3B ) 133 

DO 9 K=1 .NH21 134 

9 H( | .J.K ) = Y ( K ) 135 

10 CONTINUE 136 

DO 30 1=1 . N04 137 

DO 30 J=1 . N04 1 38 

DO 30 K = 1 . NH2 1 1 39 ' 

30 HH ( 1 • J . K ) =HI I . J.K ) 1«0 

WRITE (9) HH 141 

DO 35 1=1 » N04 1 42 

DO 35 J=1.N04 143 

DO 35 K = 1 .NH21 144 

35 HH< I .J. KI=H< I .N04+J.K1 145 

WRITE (9) HH 146 

DO 40 1=1 ,N04 147 

DO 40 J=1 . N04 148 

DO 40 K = 1 .NH21 149 

40 HH ( I • J iK ) =H< I . J .K ) 150 

WRITEI91 HH 151 

00 45 1=1 ,N04 152 

DO 45 J=1 . N04 1 53 

DO 45 K=! .NH21 1 54 

45 HH< I . J.K )=H(N04+| .J.K 1 155 

WRITE<9) HH 15 6 

DO 50 1=1 . N 04 1 57 

DO 50 J=1 « N04 158 

DO 50 K = 1 .NH21 159 

50 HH ( I . J.K 1 =H (N04+I .N04+J.K ) 160 

WRITEtS) HH 161 

DO 55 1=1 . N04 162 

DO 55 J=1 « N04 163 

DO 55 K = 1 . NH2 1 164 

55 HH ( 1 . J.K 1 =H(N04+I . J.K) 165 

WRITE (9) HH 166 

REWINO 9 167 

DO 15 K = 1 . NH2 1 168 

DO 15 1=1 .N21 1 69 

15 HN21 ( I .K )=H( I .N21 .K ) 170 

RETURN 171 

END 1 72 
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c 

c 

c 

c *********** ********************************************************* ******* 
c#*** *********************************************** ************************ 

C THE FOLLOWING IS THE GETPHI OVERLAY* WHICH COMPUTES THE POTENTIAL MESH. 

C IT REPLACES CHUNKS OF DENSITY STORED ON DISK FILES 1*2*5 AND 6 WITH 
C CORRESPONDING CHUNKS OF THE POTENTIAL MESH. IT IS CALLED ONCE PER TIME 
C STEP. 

OVERLAY (GF I LE *5*0 ) 

PROGRAM GETPHI 

C THIS OVERLAY SOLVES FOR THE POTENTIAL MESH (DIMENSIONED N/2 BY N/2 BY 

C NH/2 ) DUE TO A DENSITY MESH (DIMENSIONED N/2 BY N/2 BY NH/2 ) BY DOING A 

C PERIODIC ANALYSIS OF THE DENSITY AND THEN A PERIODIC SYNTHESIS OF THE 
C PRODUCT OF THE TRANSFORMED GREENS FUNCTION (DIMENSIONED (N/2+1) BY (N/2+1) 
C BY (NH/2+1 ) > AND THE TRANSFORMED DENSITY. FORMALLY SPEAKING* EACH OF THE 
C TRANSFORMS (EXCEPT THE COSINE ANALYSIS OF THE GREENS FUNCTION* WHICH IS 
C PERFORMED IN THE GETH OVERLAY) REQUIRES AN ARRAY DIMENSIONED N BY N JY 
C NH. TO REDUCE CORE STORAGE THIS OVERLAY PERFORMS THESE TRANSFORMS IN 
C CHUNKS BY THE ALIGNMENT OF FOUR SMALLER ARRAYS NAMED RHOl • RH02 • RH03 * AND 
C RH04* EACH OF WHICH IS DIMENSIONED N/4 BY N/4 BY NH/2. THE CHUNKS OF THE 
C LOWER HALF <1 .LE. 2 .LE. NH/2) CF THE EXTENDED ARRAY NOT IN CORE AT ANY 
C ONE TIME ARE STORED ON DISK FILES 1 THROUGH 8. THE FOLLOWING ARE TWO TOP 

C VIEWS OF THE LOWER HALF OF THE EXTENDED ARRAY. BOTH OF THESE VIEWS 

C DESIGNATE THE CHUNKS AS I ROW ANQ JCOLUMN. I ROW 1 AND 2 OF JCOLUMN 
C 1 AND 2 CONSTITUTE THE ACTIVE MESH. IN THE OIAGRAM ON THE LEFT THE 
C NUMBERS WITHIN THE CHUNKS OF JCOLUMN I AND 2 INDICATE THE DISK FILES ON 
C WHICH THOSE CHUNKS ARE STORED. (NO DISK STORAGE IS REQUIRED FOR JCOLUMN 3 
C OR 4.) REFERRING TO THE DIAGRAM ON THE RIGHT. THE NUMBERS WITHIN THE 
C CHUNKS ARE THE ORDF.R IN WHICH CHUNKS OF THE TRANSFORMED DENSITY ARE 
C MUTIPLIED (ELEMENT BY ELEMENT) BY THE APPROPRIATE PORTION OF THE 
C TRANSFORMED GREENS FUNCTION WHICH HAS BEEN READ FROM DISK FILE 9 INTO 
C ARRAY HH (N/4 *N/4 *NH/2+t ) . (AN EXCEPTION IS THE SET OF TRANSFORMED GREENS 
C FUNCTION BOUNDARY VALUES FOR I=N/2+l AND J=N/2+l WHICH REMAIN AT ALL TIMES 
C IN COMMON IN THE ARRAY HN2 ! ( N/2 ♦ 1 * NH/2+ 1 ) • ) A PLUS IN A CHUNK INDICATES 
C THAT NEW VALUES MUST' BE READ INTO ARRAY HH BEFORE THAT CHUNK IS MULTIPLIED 
C BY HH. THIS SYSTEM MINIMIZES PERIPHERAL PROCESS TIME BY UTILIZING THE 
C PERIODICITY OF THE TRANSFORMED GREENS FUNCTION. 

C 

C 

C 

C TWO TOP VIEWS OF LOWER HALF OF EXTENDED MESH ( N BY N BY NH/2) - I ROW 1 

C AND 2 OF JCOLUMN 1 AND 2 CONSTITUTE THE ACTIVE MESH (N/2 BY N/2 BY 

C NH/2). THE DIRECTIONS ARE X(t) AND OMEGAX ( I ) - DOWN ON PAGE* 

C Y ( J ) AND OMEGAY(J) - TO RIGHT ON PAGE. Z(K) AND OMEGA Z ( K ) - OUT OF 

C PAGE. 

C 
C 
C 
C 
C 
C 

c 

C I ROW = 1 

c 
c 

C I ROW=2 

C 
C 

C I ROW=3 

C 
C 

C I R0W=4 

c 
c 

C DISK FILES ON WHICH CHUNKS ORDER IN WHICH CHUNKS ARE 

C ARE STORED MULTIPLIED BY APPROPRIATE 

PORTION OF TRANSFORMED 
GREENS FUNCTION 


JCOLUMN 
12 3 4 

***************** 
* * * * * 


JCOLUMN 
12 3 4 

***************** 
* + * + * * * 

*1*5* * * 

***************** 
* * * * * 

I ROW = 1 

*1 *3*2*4* 
***************** 
* + * + * * * 

*2*6* * * 
***************** 
* * * * * 

I ROW = 2 

* 9 *11 *10 *12 * 

***************** 

* + * * * * 

*3*7* * * 

***************** 
* * * * * 

I ROW = 3 

*7*5*8*6* 
***************** 
* + * * * * 

*4*3* * * 

***************** 

I ROW = 4 

*15 *13 *16 *14 * 
***************** 


173 

174 
1 75 
1 76 
1 77 
1 78 
1 79 
1 60 
181 
182 
1 83 
1 84 

185 

186 
1 87 
1 68 
1 89 
1 90 

191 

192 

193 
1 94 
1 95 
1 96 
1 97 
1 98 
1 99 
200 
201 
202 
2C3 

204 

205 

206 

207 

208 

209 

210 
21 1 
212 

213 

214 

215 

216 
217 
21B 

219 

220 
22 1 
222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 

241 

242 

243 
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C 244 

C0MM0N/ALLC0M/N*N02*N2I * N04 . N34 ♦ NH , NH02 * NH2 1 * I 2A * 1 2B * I3A • I 3B 245 

COMMON/TRANCOM/RHOl (32*32.16) *RH02 (32 *32 • 16 ) *RH03 ( 32 * 32 * 1 6 ) * 246 

1 RH04 (32*32* 16 ) *HH( 32.32* 17) 247 

COMMON/HN21 COM/HN21 ( 65* 1 7 ) 248 

C THE INITIALIZING OVERLAY OR STAR ADVANCING OVERLAY STORES THE DENSITY 249 

C CHUNKS OF IROW 1 AND 2 FOR JCOLUMN=l ON DISK FILES 1 AND 2 RESPECTIVELY 250 
C AND FOR JCOLUMN=2 ON DISK FILES 5 AND 6 RESPECTIVELY. THE GETPHI OVERLAY 251 
C REPLACES THE DENSITY ON THESE DISK FILES WITH THE CORRESPONDING VALUES OF 252 
C POTENTIAL WHICH ARE THEN USED IN THE STAR ADVANCING OVERLAY. THIS IS 253 

C ACCOMPLISHED THROUGH CALLING SUBROUTINES ANL X ( JCOLUMN ) . ANLSYN(IROW) AnD 254 
C SYNX (JCOLUMN) AS OETAILED BELOW. 255 

C 256 

C 257 

C SUBROUTINE ANLX (JCOLUMN) READS RESPECTIVELY IROW 1 AND 2 FROM THE 258 

C FOLLOWING DISK FILES - 1 AND 2 FOR JCOLUMN* 1« - 5 AND 6 FOR JC0LUMN=2. 259 

C IT THEN PERFORMS A PERIODIC ANALYSIS IN THE X DIRECTION OVER JCOLUMN FOR 260 

C 1=1 *N AND WRITES ThE RESULTS RESPECTIVELY FOR IROW 1*2*3* AnD 4 ON THE 261 

C FOLLOWING DISK FILES - 1*2*3* AND 4 FOR JCOLUMN- 1 * - 5*6*7* AND 8 FOR 262 

C JCOLUMN=2* 263 

CALL ANLX ( 1 ) 264 

CALL ANLX (2) 265 

C SUBROUTINE ANLSYN(IROW) READS RESPECTIVELY JCOLUMN 1 AND 2 FROM THE 266 

C FOLLOWING DISK FILES - 1 ANO 5 FOR IROW* I* - 2 AND 6 FOR IR0W=2. - 3 AND 267 

C 7 FOR I ROW = 3 * - 4 AND 8 FOR IR0W = 4. IT THEN PERFORMS A PERIODIC ANALYSIS 268 

C IN THE Y DIRECTION OVER IROW FOR J=1.N. FOR EACH CHUNK It THEN PERFORMS A 269 
C PERIODIC ANALYSIS IN THE Z DIRECTION FOR K=1*NH* ELEMENT BY ELEMENT 270 

C MULTIPLICATION WITH A SIMILARLY SHAPED CHUNK OF THE TRANSFORMED GREENS 271 

C FUNCTION AND THEN A PERIODIC SYNTHESIS IN THE Z DIRECTION FOR K=1.NH. THE 272 
C RESULT FOR K=l*NH/2 IS THEN PERIODICALLY SYNTHESIZED IN THE Y DIRECTION 273 
C OVER IROW FOR J=1*N* THIS LAST RESULT FOR JCOLUMN 1 AND 2 IS WRITTEN 274 

C RESPECTIVELY ON THE FOLLOWING DISK FILES - 1 AnD 5 FOR IR0W=1* - 2 AND 6 275 

C FOR I R0W = 2 * - 3 AND 7 FOR IR0W=3* - 4 AND 8 FOR IR0W = 4. THE ORDER IN 276 

C WHICH ANLSYN IS CALLED FOR IROW I THROUGH 4 MINIMIZES READING FROM DISK 277 

C FILE 9 OF CHUNKS OF THE TRANSFORMED GREENS FUNCTION AS MENTIONED ABOVE. 27B 

CALL ANLSYN ( 1 ) 279 

CALL ANLSYN ( 3 ) 280 

CALL ANLSYN (2) 281 

CALL ANLSYN (4) 282 

C SUBROUTINE SYNX ( JCOLUMN ) READS RESPECTIVELY IROW 1*2*3* AND 4 FROM THE 283 

C FOLLOWING DISK FILES - 1.2*3* AND 4 FOR JCOLUMN= 1 • - 5*6.7, AND 8 FOR 284 

C JC0LUMN=2* IT THEN PERFORMS A PERIOOIC SYNTHESIS IN THE X DIRECTION OVER 285 
C JCOLUMN FOR J=1 *N. IT THEN WRITES THE RESULT RESPECTIVELY FOR 1 ROW • 1 AND 286 
C 2 ON THE FOLLOWING DISK FILES - 1 AND 2 FOR JCOLUMN* I' - 5 AND 6 FOR 287 

C JC0LUMN=2. 288 

CALL SYNX ( 1 ) 289 

CALL SYNX (2) 290 

RETURN 291 

END 292 

SUBROUTINE ANLX ( JCOLUMN ) 293 

C0MM0N/ALLC0M/N.N02 *N21 . NO 4 « N34 • NH • NH02 * NH2 1 * I2A* I2B* I3A* I 3B 294 

COMMON/TRANCOM/RHOI (32.32* 1 6 ) *RH02 ( 32 • 32 . 1 6 ) • RH03 ( 32*32* 16) • 295 

1 RH04 ( 32* 32 • 16 ) *HH( 32*32 • 17 ) 296 

COMMON Z ( 1 025 ) • Y ( 1 025 ) 297 

IF (JCOLUMN. EO*2 ) GO TO 2 298 

READ ( 1 ) RHO 1 299 

REWIND 1 300 

READ < 2 ) RH02 301 

REWIND 2 302 

GO TO 3 303 

2 CONTINUE 304 

READ (5 ) RHO 1 305 

REWIND 5 306 

READ ( 6 ) RH02 307 

REWIND 6 308 

3 CONTINUE 309 

CALL GETSET < 3. I2A ) 310 

DO 10 K=1 .NH02 31 1 

DO 10 J = 1 .N04 312 

00 51=1 . N04 313 

Z< I >=RH01 ( I • J «K ) 314 

Z (N04+I )=RHQ2 ( I * J*K ) ' 315 
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Z(N02+I 1=0. 316 

5 Z1N34+I )=0. 317 

CALL FTRANS<3. I2A ) 318 

00 10 1=1 .N04 31 9 

RHOl ( 1 . J.K ) =Y ( I ) 320 

RH02 < I . J.K ) =Y (N04+1 ) 321 

RH03 < I .J.K)=Y(N02+I ) 322 

10 RH04 ( I . J.K ) =Y (N34+ I ) 323 

IF ( JCOLUMN.EO.2 ) GO TO 12 324 

WRITS ( 1 ) RHOl 325 

Rewind i 326 

WRITE (2) RH02 327 

REWIND 2 328 

WRITEI3) RH03 329 

rewind 3 330 

WRITE<4> RH04 331 

REWIND 4 332 

GO TO 15 333 

12 CONTINUE 334 

WRITE <5 ) RHOl 335 

REWIND 5 336 

WRITE (6 I RH02 337 

REWIND 6 338 

WRITE (7) RH03 339 

REWIND 7 340 

WRITEIB) RH04 341 

REW I NO 8 342 

15 RETURN 343 

END 344 

SUBROUTINE ANLSYN(IROW) 345 

COMMON/ALLCOM/N . N02 . N2 1 .N04.N34.NH.NH02.NH21 . 1 2A , I2B. I 3A . I3B 346 

COMMON/TRANCOM/RHOl <32.32. 1 6 ) « RH02 < 32 « 32 . 16 ) . RH03 < 32.32. 16 ) • 347 

1 RH04 ( 32. 32 . 16 ) .HH< 32 .32 . 17 ) 348 

C0MM0N/HN21 C0M/HN21 <65. 1 7 ) 349 

COMMON Z < 1 025 ) . Y < 1 025 ) 350 

GO T0<1.2.3.4) 1 ROW 351 

1 CONTINUE 352 

READ! 1 ) RHOl 353 

REWIND 1 354 

READ <5 ) RH02 355 

REWIND 5 356 

GO TO 5 357 

2 CONTINUE 358 

READ <2 ) RHOl 359 

REWIND 2 360 

READ ( 5 ) RH02 361 

REWIND 6 362 

GO TO 5 363 

3 CONTINUE 364 

Read < 3 > rhoi 365 

REWIND 3 366 

READ <71 RH02 367 

REWIND 7 368 

GO TO 5 369 

4 CONTINUE 370 

READ < 4 ) RHOl 371 

REW INC 4 372 

READ <81 RH02 373 

REWIND 0 374 

5 CONTINUE 375 

CALL GETSET < 3 . 1 2 A ) 376 

DO 10 K=1 ,NH02 377 

DO 10 1=1 . N04 378 

DO 7 J = 1 « N04 379 

Z < J)=RH01 < I .J.K 1 380 

Z (N04 + J )=RH02 < I . J.K ) 381 

Z < N02 + J ) =0. 382 

7 Z ( N34 + J ) = 0 . 383 

CALL FTRANS < 3 * I 2a > 384 
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00 10 J = 1 . N04 385 

RH01 f I » J»K ) = Y( J 1 386 

RM02 ( I . J.K ) = Y (N04+J ) 387 

RH031 1 . J.K > =Y (N02+J > 368 

10 RH04 ( I . J.K ) = Y1N34+J ) 389 

GO TO130.49.75.75 ) I ROW 390 

49 CONTINUE 391 

30 CONTINUE 392 

R£AD(9) HH 393 

50 JC0LUMN=1 394 

DO 70 1=1 • N04 395 

00 70 J=1 . N04 396 

DO 52 K=1.NH02 397 

Z ( K ) sRHO 1 ( I t JtK ) 398 

52 Z ( NH02+K 1=0# 399 

CALL GETSET13. I3A 1 400 

CALL FTRANS < 3. I3A 1 401 

IEMR0W.NE.31 GO TO 300 402 

IEM.NE.11G0 TO 300 403 

LL=J 404 

GO TO 200 405 

54 DO 70 K=l i NH02 406 

70 RN01 ( I • J.K ) = Y <K 1 <07 

GO TO 100 408 

74 CONTINUE 409 

READ (9) HH 410 

75 JC0LUMN*2 411 

DO 95 1=1 • N04 4 1 2 

DO 95 J = 1 . N04 413 

DO 77 K = 1 . NH02 414 

Z<K)=RH02 ( 1 . J.K ) 415 

77 Z ( NH02+K ) =0 . 416 

CALL GETSET13. 1 3A ) 417 

CALL F T R A NS ( 3 . I3A ) 418 

I F ( | ROW .NE . 3 1 GO TO 300 419 

IFU.NE.llGO TO 300 420 

LL=N04+J 421 

GO TO 200 422 

79 DO 95 K = 1 « NH02 423 

95 RH02 ( 1 . J.K ) = Y (K ) 424 

GO TO 125 425 

100 JC0LUMN=3 426 

DO 120 1=1 «N04 427 

DO 120 J=1«N04 428 

DO 101 K=1.NH02 429 

Z <K I =RH03 < I . J «K I 430 

1 0 1 Z(NH02+K)=0. 431 

' CALL GE TSE T13.I3A) 432 

CALL FTRANS13. I3A 1 433 

GO T0( 1 03. 105. I 07. 1 15 1 I ROW 434 

103 1 F ( J . NE • 1 ) GO TO 300 435 

LL=! 436 

GO TO 200 437 

105 I F { j . NE . 1 1 GO TO 300 438 

LL=N04+ 1 439 

GO TO 200 440 

107 IF ( I .NE. 1 .AND. J.NE. 1 1G0 TO 300 441 

IF< | .EQ.l .AND.J.EQ. 1 >G0 TO 111 442 

1 F ( I . EQ . 1 ) GO TO 109 443 

LL= I 444 

GO TO 200 445 

109 LL=J 446 

GO TO 200 447 

1 1 1 LL =N2 1 448 

GO TO 200 449 

115 IFCJ.NE.l ) GO TO 300 450 

LL=N04+! 451 

GO TO 200 452 

117 DO 120 K=1.NH02 453 

120 RH03I I . J.K)=Y(K ) 454 

GO TO(74. 74. 400.390) I ROW 455 

125 JC0LUMN=4 456 
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DO 145 1=1 • N04 


457 


DO 145 J= 1 iN04 


458 


DO 127 K= 1 * NH02 


459 


Z(K)=RH04 ( I , J.K ) 


460 

12? 

Z (NH02 + K ) =0. 


461 


CALL GETSET ( 3 , I 3 A ) 


462 


CALL FTRANS( 3. I 3A ) 


463 


I F < I ROW • NE • 3 ) GO TO 300 


464 


IF ( I .NE. 1 )GO TO 300 


4 65 


LL=N04+ J 


466 


GO TO 200 


467 

129 

DO 145 K= 1 * NH02 


468 

1 45 

RH04 ( I.J,K)=Y(K) 


469 


GO TO (400*400,49.49) I ROW 


4 70 

200 

DO 205 K=2.NH02 


471 


Z(K)=Y(K)*HN21 ( LL * K ) 


4 72 

20 5 

Z(NH02+K)=Y(NH02+K)*HN21 (LL.K) 


4 73 


Z ( 1 )=Y( 1 )*HN21 (LL* 1 ) 


474 


ZINH21 ) - Y ( NH2 1 )*HN2l (LL.NH21 ) 


475 


GO TO 310 


4 76 

300 

DO 305 K=2.NH02 


477 


Z(K)=Y(K)*HH( I , J,K) 


478 

305 

Z(NH02 + K )=Y(NH02*fK)#HH( ! • J.K ) 


4 79 


Z( 1 ) — Y ( 1 )*HH< 1 , J, 1 ) 


480 


Z (NH21 ) = Y (NH21 )*HH< I , J.NH21 ) 


481 

310 

CALL GETSET (4 . I3A ) 


482 


CALL FTRANS ( 4 , I 3A ) 


483 


GO TO (54* 79. 1 17.129) JCOLUMN 


484 

390 

REWIND 9 


485 

400 

CALL GETSET (4 . \ 2A ) 


4 86 


DO 4 10 K = 1 * NH02 


407 


DO 410 1 = 1 * N04 


488 


DO 405 J= 1 • N04 


489 


Z( J)=RHOl ( I « J «K > 


4 9 0 


Z (N04+J ) =RH02 ( I . J.K ) 


491 


Z (N02+J >=RH03 ( I .J.K ) 


492 

405 

Z (N34+J )=RH04 ( I .J.K ) 


493 


CALL FTRANS ( 4 • I 2 A ) 


494 


DO 410 J=1 .N04 


495 


RH01 ( I . J.K)=Y< J) 


496 

410 

RH02 ( 1 . J.K )=Y(N04+J) 


497 


GO T0(415. 420. 425*430) I ROW 


498 

415 

CONTINUE 


499 


WRITE( 1 ) RHOl 


500 


REWIND 1 


501 


WRITE (5) RH02 


502 


REWIND 5 


503 


GO TO 435 


504 

420 

CONTINUE 


505 


WRITE12) RHOl 


506 


REWIND 2 


507 


WRITE (6) RH02 


500 


REWIND 6 


509 


GO TO 435 


510 

425 

CONTINUE 


511 


WRITE<3> RHOl 


512 


REWIND 3 


513 


WRITE<7) RH02 


514 


REWIND 7 


515 


GO TO 435 


516 

430 

CONTINUE 


517 


WR I TE < 4 ) RHOl 


518 


REWIND 4 


519 


WRITE <8) RH02 


520 


REWIND 8 


521 

435 

RETURN 


522 


END 


523 


SUBROUT l NE SYNX < JCOLUMN ) 


524 


COMMON/ALLCOM/N * N02 • N2 1 * NO 4 . N34 * NH.NH02.NH2 1 

, I 2A, I 2B. I3A, I3B 

525 


common/trancom/rhoi ( 32 . 32 . 16 > .RH 02 . 02 . 32. 16 > 

, RH03 (32,32,16). 

526 

1 

1 RH04 (32.32, 16 ) .HH( 32. 32 . 17 ) 


527 


COMMON Z ( 1025 ) ,Y ( 1 025 ) 


528 


IF(JC0LUMN.EQ.2 ) GO TO ! 


529 



APPENDIX 


READU) RHOI 530 

REWIND 1 531 

Read (2) RH02 532 

REWIND 2 533 

RE ADO) RH03 534 

REWIND 3 535 

Read <4) rho4 536 

REWIND 4 537 

GO TO 2 53e 

1 CONTINUE 539 

REAOI5) RH01 540 

REWIND 5 54 1 

read < 6 ) rho2 5 42 

REWIND 6 543 

Read (7) rho3 544 

REWIND 7 545 

READ ( 8 ) RH04 540 

REWIND B 54 7 

2 CONTINUE 54 8 

4 CALL GETSET14 , I2A ) s49 

DO 10 K=1 ,NH02 550 

DO I 0 J=1 . N04 55 , 

DO 5 1=1 « N04 552 

Z ( I )=RH01 I I .J.K > 553 

ZIN04+I )=RH02< I .J.K) 554 

2 IN02 + I >=RH03 ( I . J .K ) 555 

5 ZIN34+I )=RH04 ( I .J.K ) 55 6 

CALL FTRANS (4, I2A ) g 57 

DO 10 I=1.N04 55 g 

RHOI I I . J.K ) =Y C I ) 5 S9 

10 RH02 ( I « J.K ) =Y (N04+1 ) g 60 

I F < JCOL.UMN.EQ. 2 ) GO TO 12 56 , 

WRITE! 1 ) RHOI 552 

REWINO 1 553 

WRITE12) RH02 554 

REWIND 2 555 

GO TO 15 555 

12 CONTINUE 567 

WRITE! 5) RHOI 55 B 

REWIND 5 569 

WRITEI 6 ) RH02 g 70 

REWIND 6 57 , 

15 RETURN 5 7 o 

END 573 
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t = 1.75 t = 2.00 t = 2.25 



t = 2.50 t = 2.75 t = 3.00 


Figure 1 Evolution of an initially balanced, infinitesimally thin disk 
of 100 000 stars with an exponential radial density variation. The 
stars have an initial velocity dispersion given by Toomre's criterion 
(ref. 5) . Time is given in units of 2n/to 0 . 
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for the two-dimensional disk system. Note that the final density at t = 3.0 approaches 
the sum of two exponentials. 
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(b) Azimuthal velocity dispersion. 




distribution 








total potential energy for the two-dimensional disk system 
Note that the ratio of rotational to potential energy is 
approaching the value of 0.14 predicted by Ostriker and 
Peebles (ref. 8) as being required for stability. 



bar instability results in a rapid expansion in the plane 










Figure 10.- Concluded. 




























Evolution of the azimuthally averaged projected surface mass density as a function 
of radius for the three-dimensional disk system. 































Evolution of the azimuthally averaged, mean radial velocity as a function of 
radius for the three-dimensional disk system. 










<b) Azimuthal velocity dispersion 












velocity for the three-dimensional disk system. 







of inertia I (divided by their initial values) for the three 
dimensional disk system. 
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Figure 19.- Evolution of various ratios of kinetic energy to 
total potential energy for the three-dimensional disk system. 
Note that the evolutions of the energy ratios are similar to 
those shown in figure 9. 




t = 2.50 t = 2.75 t = 3.00 


Figure 20.- Evolution of an infinitesimally thin exponential disk 
with a self-consistent exponential core component. Note that 
the evolution is considerably less violent than that displayed 
in figure 1 . 
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disk-core system. Note that there is practically no change in the core structure, 
the bar instability in the disk causes an increase in the density at larger radii. 


































Comparison of the circular velocity V* and the mean azimutha] 







[ure 28.- Evolution of various ratios of kinetic energy to 
potential energy for the two-dimensional disk-core system 








■* M 


O 

m 

I! 


l r\ 


II 


a) 

p 

O 

0 

1 

to 


•H XI 
T3 P 

i— I 0) 

(a x: 
C P 

o 

•h x: 

to P 




c 

•»H 




<D 

3 




B 





•H 

T3 




■a 

<U 




i 

P 




0) 

m 




(U 

a 




P 

s 




x; 

o 




p 

o 



• 

<d 

c 

• 


c 

x: 

QJ 

o 


o 

P 

sz 

r— 


•rH 


3 


o 

4J 

<p 


<u 

• 

U 

o 

>1 


CM 

0) 


p 

a 

II 

-n 

c 

■H 

cn 

0 

o 

rH 

*H 

>H 

•H 

•H 

*W 


0* 

p 

XI 




3 

ra 

c 


N 

tH 

p 

•H 


1 

o 

to 



X 

> 

<1> 

<u 

g 


, 


■ — i 

o 


(0 

d) 

X 

X 



XI 

<0 

CO 



p 






p 

e 



ip 

(0 

<u 



o 

B 

p 




tt) 

CO 


to 
s 
<d a> 


>1 

to 


-h x: j* 
> P to 


<d >a 

4 J 

O rH 

z <o 
c 
o 

• -H 

e to 

Q) C 
p a> 
a> to e 
>1 -h 


ai 

Tl 


cn 


CT> 

(N 


P 

3 

O' 

•H 

b 


to -O 


81 



















Figure 30.- Evolution of the three-dimensional exponential disk-core 
system viewed in the equatorial (x-y) plane. Note the development 
of the comparatively weak spiral structure. 
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Figure 31 Evolution of the velocity distribution as a function of radius for 
, the three-dimensional disk-core system. 


















(c) Axial 











disk-core system. Note that there is practically 





r, kpc 

mass density for 




A Core 



Figure 35.- Evolution of the azimuthally averaged velocity dispersions as a function of 

radius for the three-dimensional disk-core system. 








Figure 35.- Continued. 








, kpc 






of inertia I (divided by their initial values) for the three' 
dimensional disk-core system. 
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t. Rotations 

Figure 38.- Evolution of various ratios of kinetic energy to 
total potential energy for the three-dimensional disk-core 


Figure 39.- Orthogonal views of the evolution of an initially non 
rotating spherical system of 1 00 000 stars (model I) . Time is 
given in units of 2 t \/& 0 . Note that the system remains spheri 
cally symmetric. 















Figure 40.- Evolution of an initially spherical system of 100 000 stars 
with an initial solid body rotation given by ft = 0.5^ o (model II). 
Note that the final state is an oblate system. 















Figure 41.- Evolution of an initially spherical system of 100 000 stars 
with an initial solid body rotation given by = 0.707^ o (model III) 
The system remains axisymmetric and acquires an oblate shape. 
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(a) Evolution with each projection containing 100 000 stars. 


Figure 42.- Evolution of initially spherical system with an initial 
solid body rotation given by fi = 0.866fl o (model IV). The 
initial energy in random motion is 1/2 that of the systems shown 
in figures 39 through 41 . Note that the system quickly forms a 
bar and has reached an essentially steady state at t = 5. 

























































(a) Evolution with each projection containing 100 000 stars. 


Figure 43.- Evolution of an initially spherical system with an initial 
solid body rotation given by = 1.159^, (model V). The initial 
energy in random motion is 1 /5 that of the systems shown in 
figures 39 through 41 . Again the system forms a bar. 
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Model II 
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t, Rotations t. Rotations 

Figure 44.- Evolution of the moment of inertia I and the angular momentum 
(divided by their initial values) for models I, II, III, and V. 




of inertia I and the an< 
tial values) for model IV 
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tential energy for model IV. 



Model I f Model II 



(3)N 



(3)N 


115 


the ordinate is the number of stars in an energy interval of 
width 0.05GMg/R. 





Number of stars per interval 5 r 5 
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Figure 51 Continued 




Figure 51 Continued. 






Figure 51 Concluded 





, MQ/kpc^ /u , MQ/kpc 
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Figure 54.- Final volume mass density variation in the axial 
direction at three cylindrical radii for models I, II, II 


a r - 5 kp 
♦ r = 10 k 



Evolution of the volume mass density variation in the axial direction 
at three cylindrical radii for model IV. 







< V 



Figure 56.- Comparison of the final circular velocity r u> 0 = \j K r (z=0) r 
with the mean azimuthal velocity <V ( j ) > for models I, II, III, 
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solid body rotation for the bar-forming model V. 




130 


radius r for model IV. Again this bar-forming model shows only a small central 
region with solid body rotation. 
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radius r for model IV. Note the nearly equal values for different colatitudes. 






0 4 8 12 16 20 0 4 8 12 IS 20 

r, kpc r, kpc 

Figure 62.- Final azimuthal velocity dispersion as a function of the 
spherical radius r for models I, II, III, and V. 





igure 64.- Final axial velocity dispersion as a function of the 
spherical radius r for models I, II, ill, and V. 







Figure 65.- Evolution of the axial velocity dispersion as a function of the 

spherical radius r for model IV. 
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Figure 67.- Projected isodensity contours for various viewing directions of model V 
at t = 3. I Q represents the peak in the projected number density of stars per 
kpc^. The contour levels have the same meaning as in the preceding figure. 
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s = 


141 


numbers of the disk files on which those chunks are stored 
and 4 do not require disk file storage. 
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